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ABSTRACT 

The  report  studies  the  constitutive  behaviour  of  highly  porous  materials.  This  study  is  driven 
by  efforts  for  the  development  and  validation  of  a  multi-phase  modelling  capability  in  DSTO, 
which  aims  at  an  enhanced  evaluation  of  blast  mitigation  by  porous  materials.  A  two-phase 
model  is  employed  for  the  description  of  a  porous  material,  for  which  experimental  shock 
velocity  data  are  widely  available  in  literature.  The  method  deriving  Hugoniots  from  the 
shock  velocity  data  is  critically  approached.  The  importance  of  the  non-equilibrium  behaviour 
analysis  is  demonstrated,  and  it  is  shown  that  the  stress  or  velocity  profiles  need  to  be 
considered,  using  relevant  experiments  and  constitutive  modelling.  It  is  argued  that  the 
apparent  anomalous  behaviour  of  the  Hugoniots  of  a  highly  porous  material  might  be  caused 
by  misinterpretation  of  the  experimental  data. 
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Constitutive  Modelling  of  the  Shock  Behaviour  of  a 

Highly  Porous  Material 


Executive  Summary 

Porous  materials  have  been  used  for  blast  and  shock  mitigation  for  a  long  time.  A  variety 
of  porous  materials  such  as  sand,  soil,  and  concrete  are  employed  as  countermeasures  to 
Improvised  Explosive  Devices  (IED)  and  in  a  range  of  activities  related  to  mitigation  of  the 
blast  and  shock  due  to  the  explosion  of  IEDs,  landmines,  submerged  mines  and  munitions 
in  storage. 

Accurate  characterisation  of  materials  is  important  for  establishing  modelling  capabilities. 
Some  porous  materials  have  been  characterised  extensively  in  numerous  laboratory  and 
underground  nuclear  tests,  mainly  in  the  USA  and  the  USSR.  The  test  results  were  used 
for  the  characterisation  of  a  wide  range  of  porous  materials  that  might  be  used  under 
extreme  pressures  in  explosive  devices  and  in  the  high  velocity  impact  environment.  Such 
characterisation  was  mainly  employed  for  the  determination  of  Hugoniots  (curves 
connecting  the  material  states  in  front  and  behind  the  shock  front);  in  turn,  the  Hugoniots 
were  used  to  derive  enhanced  equations  of  state.  The  Hugoniots  were  obtained  with  the 
impedance  mismatch  method.  Hugoniots  are  sometimes  the  only  data  available  for 
derivation  of  the  equations  of  state  that  allow  hydrocode  modelling  and  evaluation  of  the 
material  state  at  high  pressures.  The  data  used  for  the  derivation  of  the  Hugoniots  are 
mainly  the  shock  velocities  in  a  tabulated  (standard)  material  and  in  a  material  under 
investigation.  Porous  materials  at  extreme  pressures  exhibit  fairly  conventional  behaviour 
of  the  Hugoniots.  However,  at  moderate  pressures  of  the  order  of  several  GPa  the  derived 
Hugoniots  for  highly  porous  materials  show  an  anomalous  behaviour  resulting  in  a 
pressure  increase  with  a  decrease  of  density.  A  physical  explanation  was  suggested  for 
this  abnormality,  which  necessitated  temperature  sensitive  effects  during  compression  of 
the  material.  However,  experimental  observations  of  shock  profiles  in  porous  materials 
demonstrate  that  the  shock  deformation  may  occur  in  several  stages,  thus,  reducing  the 
thermal  contribution  into  the  energy  absorbed  by  the  material.  Although  the  thermal 
contribution  might  be  overestimated,  the  nature  of  the  apparent  abnormality  of  the 
Hugoniot  has  not  yet  been  completely  explained. 

The  present  report  describes  the  behaviour  of  a  highly  porous  material  at  shock 
compression  using  a  constitutive  two-phase  model.  Conventional  equations  of  state  are 
employed  for  the  material  components,  which  have  no  abnormalities.  The  present 
modelling  enables  a  good  description  of  available  shock  velocity  experimental  data  for  a 
highly  porous  copper.  The  present  work  has  shown  that  experimental  data  on  stress 
profiles  and  constitutive  analysis  of  the  behaviour  of  porous  materials  are  necessary  for 
the  material  characterisation  and  evaluation  of  the  material  response  to  blast  and  shock. 
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1.  Introduction 


The  characterisation  of  porous  materials  has  attracted  the  attention  of  researchers  because  of 
the  possibility  of  expanding  the  range  of  state  parameters  into  the  high-temperature  region. 
Testing  at  extreme  pressures  was  possible  during  the  extensive  underground  nuclear  testing 
programs  in  the  United  States  and  the  former  Soviet  Union  in  the  1950s.  In  particular, 
researchers  from  the  Russian  Federal  Nuclear  Centre  (RFNC)  of  Arzamas,  (former  Sarov, 
Nizhnij  Novgorod  region)  reached  the  Terapascal  range  of  pressures  using  test  assemblies 
embedded  in  the  rock  surrounding  a  detonating  nuclear  device  [1].  The  researchers  tested 
various  solid  and  porous  metals  and  organic  compounds.  Some  of  the  tests  were  also  carried 
out  in  laboratory  conditions  in  the  Gigapascal  range  of  pressures,  using  high  explosive  (HE) 
planar  and  spherical  wave  generators.  To  achieve  simplicity  in  experimental  set-up  the 
researchers  used  the  Impedance  Mismatch  (IM)  method  (see  [2];  this  method  is  named  the 
reflection  method  [3]  in  the  Russian  literature).  The  tests  allowed  researchers  to  measure 
shock  wave  velocities  in  a  standard  and  in  a  material  under  investigation  simultaneously. 
When  having  IM-processed  the  data,  the  researchers  inferred  pressure  and  density  in  a  tested 
material  behind  a  given  shock  wave.  This  enabled  them  to  draw  both  a  Hugoniot  of  the  tested 
material  in  the  pressure-density  space  and  a  line  in  the  shock  velocity-particle  velocity 
coordinates,  which  is  frequently  also  called  a  Hugoniot.  Such  the  data  have  been  derived  for  a 
large  number  of  metals  and  porous  substances  (see,  for  example,  [2,  4,  5]).  Data  for  porous 
copper,  which  are  of  specific  interest  to  the  present  work,  have  been  reported  in  [6-21].  The 
present  work  focuses  on  an  analysis  of  data  for  porous  copper  with  an  initial  density  poo  that  is 
approximately  a  quarter  of  the  solid  copper  density  po  (the  porosity  rn  =  po/poo  =  4).  Data 
relevant  to  this  porosity  were  reported  in  [6-12].  Only  a  few  data  points  for  this  material  have 
been  found  in  US  literature  sources  [11-12] .  Thus,  the  RFNC  Russian  team  has  published  the 
majority  of  the  data,  and  the  majority  of  the  experimental  points  obtained  in  laboratory 
studies  have  been  reported  in  paper  [6],  This  data  attracts  significant  attention  from  modellers 
because,  i)  it  is  quite  typical  of  this  class  of  data  for  the  family  of  porous  materials;  and  ii)  it  is 
a  sound  representative,  the  IM-processing  of  which  results  in  the  anomalous  feature  of  the 
Hugoniot  (increase  of  pressure  with  a  decrease  of  density).  These  Hugoniots  were  used  as 
reference  curves  for  derivation  of  a  large  variety  of  Equations  Of  State  (EOSs),  which  invoke 
rather  complex  high  temperature  excitation  mechanisms  in  order  to  describe  the  abnormality 
[8, 9, 21]  occurring  in  the  pressure  range  below  50  GPa.  It  should  be  mentioned,  however,  that 
a  paper  [22]  points  out  that  the  necessity  of  considering  some  high-temperature  mechanisms 
might  be  overestimated,  when  designing  EOSs  in  the  pressure  range  below  hundreds  of  GPa. 
The  test  set-up  used  in  publications  [6-12]  enabled  the  researchers  to  draw  so-called  (D-D)- 
dependencies  that  link  the  shock  velocity  in  a  tested  material  with  the  shock  velocity  in  a 
standard  material  at  the  same  loading  conditions.  Several  examples  of  such  (D-D)-graphs 
were  drawn  in  [7].  Unfortunately,  these  dependencies  were  not  detailed  in  a  table  form,  so  an 
accurate  (D-D)-data  cannot  be  recovered  from  the  paper.  Even  in  publications  where  some 
isolated  (D-D)-points  can  be  found  (e.g.,  [17])  the  authors  report  on  a  reduced  (IM-processed) 
data  using  'a  calculated  law  of  attenuation  of  shock  waves',  which  does  not  allow  the 
experimental  conditions  and  the  test  records  to  be  recovered  with  confidence.  From  the  (D-D)- 
dependencies,  the  (P-p)-  and  (D-U)-Hugoniots  (here  P  is  pressure,  p  -  density,  and  U  - 
velocity  behind  the  shock  front)  have  been  derived  with  the  IM-technique.  The  former  (P-p)- 
dependence  is  the  Hugoniot  actually  sought  and  the  (D-U)-line  is  usually  called  a  Hugoniot  as 
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well.  Using  these  dependencies,  EOSs  are  derived  by  invoking  various  physical  mechanisms 
of  the  material  response  to  compression.  One  serious  simplification,  when  using  the  IM- 
technique,  is  the  jump  approximation  for  a  shock  wave  in  porous  material.  This 
approximation  is  accurate  enough  in  solid  metals  loaded  by  strong  shock  waves;  however, 
even  metals  subject  to  moderate  shock  waves  exhibit  non-equilibrium  behaviour.  Porous 
materials  exhibit  kinetic  response  that  is  caused  both  by  a  compression  precursor  and  by  non¬ 
equilibrium  (kinetic)  behaviour  of  the  main  shock  wave.  Therefore,  a  complete  temporal 
record  of  stress  or  velocity  could  give  proper  appreciation  of  the  state  achieved  during  shock 
loading.  However,  such  experimental  information  is  not  readily  available  and  is  reported  for 
quite  a  limited  range  of  porous  materials  and  for  limited  porosity  and  load  ranges  (see,  for 
example,  [14-16,  23-24]). 

The  publications  of  the  RFNC  scientists  on  the  shock  compression  experiments  do  not  provide 
detailed  information,  such  as  an  exact  set-up  description,  parameters  of  the  test  assemblies 
and  unprocessed  recorded  shock  velocities  for  standard  and  tested  materials.  Instead,  the 
publications  deliver  post-processed  Hugoniots  in  (P-p)-  and  (D-U)-coordinates,  which  have 
somehow  been  corrected  in  order  to  take  into  account  various  factors  such  as  the  shock  wave 
attenuation  during  its  propagation  in  the  experimental  set-up  (for  example,  the  attenuation 
correction  law  mentioned  above),  non-planarity  of  impact,  curvature  of  the  shock  front,  etc. 
The  lack  of  specifications  in  these  set-up  descriptions  and  in  the  reported  results  was  noted  by 
Nellis  et  al  in  [22].  However,  even  recent  publications  [1,  4,  25-27]  aimed  at  an  extensive 
description  of  the  test  conditions  and  set-ups  of  the  laboratory  experiments  from  [6-12]  also  do 
not  give  complete  details;  this  issue  will  be  discussed  in  Section  2. 

The  present  work  employs  a  two-phase  constitutive  model  described  elsewhere  [24]  and 
briefly  outlined  in  Section  3.  This  model  utilizes  conventional  equations  of  state  for  gaseous 
and  solid  phases  of  a  porous  material.  The  present  work  attempts  to  estimate  experimental 
set-ups  as  closely  as  possible  to  original  ones  from  the  available  published  information  for  the 
tests  on  the  porous  (m  =  4)  copper  and  to  simulate  the  data  [6-12]  using  the  constitutive  model 
[24].  The  numerical  results  describe  the  experimental  velocity  data  (the  (D-D)-  and  (D-U)- 
dependencies)  quite  well.  In  doing  so,  however,  the  Hugoniots  in  (P-p)-coordinates  obtained 
from  the  simulated  (D-D)-results  by  the  IM-processing  exhibit  no  abnormalities.  It  has  been 
shown  that  a  possible  reason  for  the  apparent  Hugoniot  abnormality  is  likely  to  be  the  use  of 
the  jump-like  (equilibrium)  approximation  for  the  analysis  of  the  shock  compression  of  the 
porous  material.  Nevertheless,  a  kinetic  description  allowing  for  an  abnormality  may  also  be 
possible.  In  that  case  the  abnormality  is  achieved,  when  modelling,  via  suppression  of  heat 
exchange  from  the  solid  phase;  this  requires  a  switch  of  the  heat  exchange  mechanisms 
resulting  in  a  pronounced  constitutive  behaviour  and  in  a  significant  non-equilibrium  behind 
the  shock  front.  Thus,  the  present  approach  demonstrates  that  the  constitutive  behaviour  of 
shock  waves  in  a  porous  material  could  be  an  essential  factor  contributing  to  the  data 
interpretation.  Therefore,  experimental  data  tracing  complete  stress  or  velocity  profiles,  such 
as  those  from  [14-16,  23-24]  would  be  of  great  assistance  in  the  evaluation  of  material  shock 
response.  Earlier,  attention  was  also  drawn  in  [28]  to  the  importance  of  data  illustrating 
constitutive  behaviour  of  porous  substances. 
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2.  Introduction  to  the  experimental  method 

The  Impedance  Mismatch  (IM)  method  is  based  on  derivation  of  a  Hugoniot  from  the  mass, 
momentum  and  energy  conservation  laws  balanced  across  the  shock  jump: 

p  (D  -  U)  =  p0  D  , 

P  =  po  D  U  ,  (1) 

E-Eo  =  1/2P-(l/p0-l/p)  , 

where,  using  an  Equation  Of  State  (EOS)  in  the  caloric  form  E=E(p,P),  the  Hugoniot  or 
Rankine-Hugoniot  curve  Ph  =  Ph(p)  can  be  obtained  from  (1).  By  an  inverse  procedure  the 
EOS  can  be  found  as  soon  as  the  Hugoniot  is  known.  In  (1)  the  initial  pressure  and  velocity  are 
assumed  to  be  negligible,  P,  p,  U,  E  are  pressure,  density,  velocity  and  internal  energy  behind 
the  shock  front,  Eo  and  po  are  internal  energy  and  density  in  front  of  the  shock  wave. 


Fig.  1  Schematic  of  the  experimental  set-up  employed  for  evaluation  of  the  material  shock  response 

A  typical  shock  response  experiment  records  the  shock  wave  velocity  D.  The  simplest  set-up 
comprises  a  screen  made  from  a  standard  material  that  has  well-known  properties  and  a 
sample  of  a  material  under  investigation.  Referring  to  Figure  1  for  the  set-up,  a  flyer  plate  (F) 
with  a  velocity  of  2Uo  made  from  a  standard  material  (S)  impacts  the  target  assembly 
consisting  a  screen  (A)  made  from  the  same  material.  The  screen  separates  the  flyer  plate  from 
samples  ( B )  and  (C)  made  from  the  material  S  and  a  tested  material  T,  which  are  positioned  in 
parallel.  The  samples  are  equipped  with  time  of  arrival  gauges  (G).  In  Fig.  1,  one  of  possible 
configurations  of  'electro-contact'  gauges  used  in  [6-10]  is  expanded,  showing  the  details 
reported  in  [21, 26]  (a  similar  description  can  be  found  in  [19];  the  gap  may  approach  1  mm  as 
described  in  [7]).  In  these  publications,  the  accuracy  of  the  technique  for  measurement  of  the 
shock  wave  velocity  in  solid  materials  is  claimed  to  be  within  1-1.5  % .  A  detailed  description  of 
the  pin  technique  has  been  reported  in  [29, 30],  including  a  discussion  on  potential  issues  and 
the  sources  of  possible  errors.  In  particular,  the  paper  [29]  is  more  cautious  about  the 
technique  accuracy  and  it  reports  data  with  accuracy  not  worse  than  5%  for  shock  velocity 
measurements  in  condensed  materials. 
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The  active  element  of  the  gauge  (pin  P )  in  Fig.  1  is  protected  by  a  cup  E  against  air  shock  that 
emerges  when  free  surface  of  the  screen  A  or  samples  B  and  C  moves,  following  the  shock 
wave  reflection  from  the  free  surface.  When  the  top  wall  section  of  the  cup  E  (the  protector 
that  will  be  referred  to  below  as  E  as  well)  has  made  contact  with  the  pin  it  transmits  a  signal 
to  an  oscilloscope.  We  denote  thicknesses  of  the  plates  F  and  A  to  be  Wf  and  Wa,  and  thickness 
of  the  samples  B  and  C  to  be  Ws.  For  the  gauge  G  we  denote  thickness  of  the  protector  E  to  be 
Wb  (0.3  mm  in  Fig.  1)  and  the  gap  width  to  be  Wg  (0. 2-0.4  mm  in  Fig.  1).  An  experiment  with 
this  set-up  determines  the  shock  velocities  Ds  and  DP  in  the  standard  and  tested  materials  at 
approximately  the  same  conditions  of  loading  (when  the  shock  wave  starts  propagating  into 
the  samples  B  and  C),  since  both  samples  are  placed  on  the  screen  A  (made  from  the  standard 
material  S). 

Below  we  illustrate  how  pressure  Po  and  material  velocity  Uo  behind  the  shock  front  in  the 
screen  can  be  calculated  as  an  intersection  point  of  a  Hugoniot  and  the  Rayleigh  line, 
providing  the  Flugoniot  of  the  standard  material  is  known  and  the  shock  velocity  Ds  in  this 
material  (sample  B)  has  been  measured. 


Fig.  2.  The  shock  propagation  (x-t)-diagram  for  the  experimental  set-up  (a),  and  the  shock  adiabat 
diagram  of  the  shock  interaction  process  (b). 

Close  examination  of  the  process  of  interaction  of  shock  waves  with  the  samples  shown  in 
Fig.  2,  shows  that  after  impact  of  the  flyer  plate  F  against  the  screen  A,  two  shock  waves  with 
the  same  velocity  Ds  go  left  and  right  into  F  and  A,  originating  from  the  contact  surface 
(. x  =  Wf).  Following  up  the  collision  event  and  propagation  of  the  right-going  wave  through 
the  screen,  the  wave  comes  first  to  the  screen-sample  interface  (A-C)  and  then  it  breaks  into 
the  shock  wave  passing  down  through  the  sample  C  and  the  rarefaction  wave  (if  impedance 
of  the  material  T  is  less  than  impedance  of  the  standard  material  S)  reflecting  up  to  the  screen 
A.  An  (x-f)-diagram  in  Fig.  2(a)  shows  the  development  of  the  shock  wave  process  in  the  test 
assembly,  where  x  is  directed  down  and  the  coordinate  origin  is  aligned  with  the  position  of 
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the  rear  surface  of  the  flyer  plate  at  the  collision  moment.  After  the  shock  wave,  travelling 
with  the  velocity  DP  through  the  sample,  has  approached  the  sample-gauge  protector  interface 
(at  x  =  Wf  +  Wa  +  Ws)  and  propagated  through  the  protector  plate  E,  the  wave  reflected  from 
the  free  surface  of  E  at  the  time  moment  t  =  h,  accelerates  the  protector  E,  so  the  protector-pin 
gap  is  closed  up  at  the  time  moment  t  =  te,  shorting  the  pin,  which  has  initially  been  located  at 
X  =  Wf  +  Wa  +  Ws  +  Wg  +  Wb. 

The  shock  interaction  process  can  also  be  tracked  using  the  shock  adiabat  diagram  drawn  in 
the  (P-U)-space  (Fig.  2(b)).  The  shock  adiabats  are  used  in  this  coordinate  system  because  a 
state  after  wave  interaction  is  more  conveniently  described  by  a  single  point  in  the  (P-U)- 
space.  The  right-going  shock  wave  propagating  into  A  gives  a  locus  of  states  (a  curve  Hs) 
connecting  the  initial  state  of  the  material  S  (state  '[0]'  at  the  origin  corresponding  to  negligible 
pressure  in  the  resting  material)  with  all  possible  states  behind  the  shock  front;  one  of  them  is 
a  state '  [1]'  at  the  impact  surface  when  colliding  the  flyer  plate  with  the  screen.  The  left-going 
shock  wave  propagating  into  F  sets  up  a  curve  that  connects  the  initial  state  of  the  flyer  plate 
'[O’]'  (it  is  under  negligible  pressure  and  its  impact  velocity  is  2Uo)  with  all  possible  states 
behind  the  shock  front;  this  curve  Hs  is  the  mirror  image  of  Hs.  Quite  a  strong  assumption  is 
made  about  the  uniformity  of  the  state  '[O']',  keeping  in  mind  that  the  flyer  plate  is  accelerated 
by  an  explosive  charge  within  a  small  stand-off  distance. 

If  the  velocity  of  the  flyer  plate  and  properties  of  the  standard  material  have  been  determined 
in  the  experiment,  then  the  curves  Hs  and  Hs  can  be  drawn  and  a  common  pressure-velocity 
state  (Po,  Uo)  (state  '[1]')  behind  the  shock  waves  with  the  velocity  magnitude  Ds  originating 
from  the  contact  (A-C)  interface  (at  x  =  Wf  in  Fig.  2(a) )  can  be  found  as  an  intersection  point  of 
Hs  and  Hs  (Fig.  2(b)).  Fiowever,  in  the  majority  of  the  experiments  the  flyer  plate  velocity  2Uo 
is  not  known  a  priori.  Instead,  the  shock  velocity  Ds  in  the  screen  material  is  measured;  more 
precisely,  it  is  measured  in  the  sample  B  attached  to  the  screen  at  x  =  Wf,  therefore,  the 
assumption  that  the  wave  is  stationary  when  propagating  within  the  assembly  A  +  B  is 
essential.  If  Ds  is  known,  then  pressure  and  velocity  behind  the  shock  front  in  S  are  connected 
by  the  momentum  equation  of  the  system  (1)  resulting  in  a  Rayleigh  line  P  =  psDsU  on  the 
diagram  in  Fig.  2(b);  here  psis  the  initial  density  of  the  standard  material.  In  this  case  the  state 
'[I]'  can  be  obtained  as  an  intersection  point  of  the  shock  adiabat  Hs  and  the  Rayleigh  line  as 
shown  in  Fig.  2(b).  Pressure  and  velocity  (Po,  Uo)  calculated  in  such  a  way  characterise  the 
shock  wave  propagating  through  the  screen  A  into  the  material  T  of  the  sample  C  under 
investigation  (a  part  of  the  A  +  C  assembly  in  Fig.  1). 

After  the  moment  t  =  h  (Fig.  2(a))  when  the  shock  wave  enters  the  sample  C,  the  shock  wave 
breaks  down  on  the  A-C  interface  (at  x  =  Wf  +  Wa  in  Fig.  2(a))  into  the  shock  wave  with 
velocity  DP  going  right  into  the  sample  C  and  the  rarefaction  wave  going  left  back  into  the 
screen  A.  The  reflection  method  (IM-method)  assumes  that  the  isentrope  (curve  Rs  in  Fig.  2(b)), 
linking  states  of  the  standard  material  with  states  behind  the  rarefaction  wave,  is  very  close  to 
the  reflected  shock  adiabat  Hs  (the  mirror  image  of  the  incident  shock  adiabat  Hs).  Similarly 
to  the  previous  consideration  for  the  standard  material  S,  states  for  the  material  T  under 
investigation  behind  the  right-going  shock  wave  propagating  with  velocity  DP  lie  on  a 
Rayleigh  line  P  =  pooDpU  in  the  (P-U)-space;  here  poo  is  initial  density  of  the  material  T.  Thus, 
state  '[2]'  behind  the  shock  front  in  the  material  T  can  be  calculated  as  an  intersection  point  of 
the  Rayleigh  line  and  the  isentrope  Rs  (approximated  by  the  curve  Hs  in  the  reflection 
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method).  This  state  (P,  U)  in  the  pressure- velocity  space  gives  a  point  of  the  shock  adiabat  Hp 
(Fig.  2(b))  for  the  material  T.  When  P  and  U  are  found,  the  density  behind  the  shock  front  is 
calculated  from  the  mass  conservation  law  (the  first  equation  of  the  system  (1)),  which 
completes  the  process  of  calculating  a  point  of  the  Hugoniot  Hp  sought.  The  shock  velocity  DP 
could  normally  be  calculated  as  Ws/  (h  -  h)  that  gives  the  average  shock  velocity,  when 
ignoring  possible  shock  wave  attenuation  (it  starts  att=h  in  Fig.  2(a)).  In  real  experiments, 
the  pin  shorting  technique  is  used  for  the  time  of  arrival  registration.  In  the  schematic  of 
Fig.  2(a),  the  shorting  of  the  pin  occurs  at  t  =  te  instead  of  at  t  =  U,  as  described  above.  Similar 
time  correction  for  the  time  of  arrival  of  shock  wave  to  the  screen-sample  interface  should  be 
made  by  taking  account  of  the  gauge  set-up  positioned  at  the  screen-sample  A-C  interface  at  x 
=  Wf  +  Wa  (a  schematic  of  the  process  for  that  gauge  is  not  shown  in  Fig.  2(a)). 

The  most  popular  issue  in  the  Flugoniot  evaluation,  which  attracts  the  attention  of  many 
researchers  (see,  e.g.,  [31,  32]),  is  an  approximation  of  the  isentrope  Rs  by  the  reflected 
Hugoniot  Hs  .  Less  attention  is  paid  to  the  non-steadiness  of  shock  waves  (kinetic  behaviour 
of  the  waves)  and  conditions  behind  the  shock  front  in  porous  materials.  Recent  advanced 
models  consider  the  constitutive  behaviour  of  porous  material;  for  example,  the  P-lambda 
model  by  Grady  [33]  implemented  in  the  CTH  hydrocode  and  used  for  analysis  of  the  shock 
wave  structure  in  porous  ceramic  powders  [34],  a  model  by  the  author  [24]  used  for  analysis 
of  shock  profiles  in  porous  aluminium  [24]  and  dry  and  water-saturated  sands  [35].  Using 
advanced  models,  it  is  possible  to  describe  the  shock  compression  of  a  porous  material  with 
more  details  than  with  the  traditional  jump-like  approximation  of  the  shock  front.  In  the 
present  work  we  will  ignore  the  issue  of  isentrope  approximation,  keeping  in  mind  that  the 
IM-reduced  data  for  the  porous  copper  reportedly  were  corrected  for  the  approximated 
mirror  adiabat  to  the  isentrope.  Thus,  the  present  study  will  mainly  concentrate  on  the  issue  of 
non-equilibrium  behind  the  shock  front.  It  should  be  mentioned  that  in  the  present  work  the 
rarefaction  waves  are  computed  directly  via  constitutive  modelling.  However,  for  analytical 
assessments  (recovery  of  the  shock  velocity  data  from  particle  velocity),  the  mirror 
approximation  will  be  used  directly,  because  the  source  publications,  reporting  the  IM- 
reduced  data,  do  not  reveal  the  nature  and  amount  of  the  corrections. 


3.  Constitutive  model 


The  constitutive  model  employed  in  the  present  work  has  been  described  elsewhere  [24]. 
Briefly,  for  a  porous  material  the  model  assumes  two  phases,  each  of  which  is  characterised  by 
density  p,  specific  entropy  s,  pressure  p,  and  temperature  T;  the  latter  two  being  dependent 
thermodynamic  variables.  Velocity  u  is  the  same  for  both  phases.  For  a  porous  material,  the 
two  phases  are  a  gas  (the  air)  and  a  solid  phase  with  negligible  strength.  Variables  for  the 
corresponding  phases  1  (gas)  or  2  (solid)  are  referred  to  by  superscripts  in  parentheses.  Both 
phases  are  assumed  to  be  compressible  with  individual  mechanical  properties.  Using  M  for 
denotation  of  mass  occupying  specific  volume  V,  the  total  mass  and  volume  are  M  =  MW  +M(2> 
and  V  =  VW  +  V(2).  Then  the  averaged  density  p  =  M/V  relates  to  corresponding  phase 
densities  pW  =  MW/  VO)  and  p(2)  =  M(2)/  V (2)  via  the  volume  and  mass  concentrations  Q  and  c  of 
the  first  (gaseous)  phase  as  follows 
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e=W)/V,  c  =  MW/M ,  p  =  0p«  +  (l-0)pP) 
Specific  internal  energy  e  and  entropy  s  are 

e  =  c  ed)  +  (1  -  c)  e(2)  ,  s  =  c  sd)  +  (1  -  c)  s(2)  , 
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(2) 

because  of  their  extensivity.  The  entropy  concentration  for  the  first  phase  j  is  introduced  in 
[24]  similar  to  the  introduction  of  the  mass  concentration.  Following  the  definitions,  we  can 
express  the  independent  phase  state  variables  via  the  averaged  ones: 

p(!)  =pc/e  ,  p( 2)  =  p(l  -  c)  /(I  -  0  )  ,  sd)  =  sx/c  ,  s(2)  =  S  (1  )/(l  ~  C)  •  (3) 

Summarising,  the  independent  variables  that  fully  describe  the  two-phase  material  are  p,  s,  u, 
c,  0,  and  j.  The  first  three  variables  are  typical  ones  involved  in  the  conservation  laws: 

dt  dxj 

8pu,  {  d(puiUj  +  pSy)  _  Q 

at  dxj  ’  1  ; 

|  d[/7^£  +  p»y]_o 

5v/ 


where  E  =  e  +  Uj-Uj/2  and  5,;  are  components  of  the  unit  tensor.  When  the  specific  internal 
energy  e  is  given  as  a  function  of  the  parameters  e  =  e(p,  s,  c,  9,  j),  the  thermodynamic  identity 
gives  p  =  p2ep  and  T  =  es. 

We  assume  that  the  internal  energy  dependencies  for  each  phase  are  given  in  the  form 

e(i)  =  e(i)(p(i),s(i))  ,  g(2)  =  e(2)(p(2),s(2))  .  (5) 

Thermodynamic  identities  for  the  phases  can  be  used  for  calculation  of  averaged  dependent 
variables  from  (2)  as  follows: 

p  =  p2ep=p2(cem  +  (l-c)el2))p=0pm  +  (l-0)p^\ 

T  =  es  =(ce(11  +  (l-c)e<2))s  =  %Ta)  +  {\- Z)T(2) .  ^ 


Considering  partial  densities  pd) p  =  M0)/V  =  p-c  and  p(2) p  =  M(2)/  V  =  p-(l-c),  we  can  associate 
the  mass  change  of  the  first  phase  with  the  rate  proportional  to  a  constitutive  function  <p.  Thus, 
the  continuity  equation  for  pd)p  gives: 


dpc  dpcUj 
dt  dxj 


-pep. 


(7) 
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Similarly,  using  the  entropy  balance  equation,  we  can  derive  the  following  balance  equation 
for  parameter  j: 


dpX  |  dPXuj 

dt  dXj 


-pa. 


(8) 


where  function  co  is  a  constitutive  function  responsible  for  the  energy  exchange  between  the 
phases.  The  constitutive  equation  for  the  volume  concentration  is  considered  in  the  following 
well-known  form  (see,  for  instance,  references  in  [24]): 


cpO  dpOiij 
dt  dXj 


-PV  ■ 


(9) 


Thus,  if  for  each  phase  the  internal  energy  functions  are  taken  in  the  form  (5),  the  specific 
energy  e(p,  s,  c,  9,  j)  can  be  calculated  from  (2)  with  the  use  of  (3).  Pressure  and  temperature 
can  be  calculated  from  (6).  The  complete  system  of  equations  for  the  model  contains  equations 
(4),  (7),  (8),  and  (9).  The  system  is  fully  defined  when  the  equations  of  state  for  each  phase  (5) 
are  pre-selected  and  the  functions  cp,  ip,  and  co  are  given.  The  entropy  increase  condition  from 
the  energy  balance  equation  imposes  the  following  restrictions  [24]  for  the  functions  cp,  ip,  and 
co  that  have  to  be  chosen  in  the  form 


cp  =  ec-cp0  ,  ip  =  ee-ip0  ,  co  =  ex-co0  ,  (10) 

where  functions  cpo,  ipo,  and  co o  are  arbitrary  non-negative  functions.  The  expressions  ec,  c’e,  and 
ex  can  be  calculated  from  (2,  5)  and  are  associated  with  the  affinities  of  the  Gibbs  potential, 
pressure  and  temperature  [24], 

Thus,  the  system  (4)  and  the  constitutive  equations  (7)-(9)  generate  a  system  describing  the 
behaviour  of  two-phase  porous  material.  The  conservation  laws  are  closed  with  the  equations 
of  state  (5)  for  each  phase,  where  parameters  of  the  phases  are  connected  with  the  averaged 
parameters  by  relations  (2)  and  (3).  The  constitutive  equations  are  closed  with  selection  of  the 
constitutive  rates  (the  arbitrary  non-negative  functions  cpo,  ipo,  and  coo)  that  will  be  specified  in 
next  Section. 


4.  Constitutive  relations  for  porous  copper 

We  choose  constitutive  functions  cpo,  ipo,  and  (0»  in  the  following  form  that  has  been  used  for 
the  constitutive  equations  for  sand  and  porous  aluminium  [24] : 

cpo  =  c(l-  c)-Ao  ,  ipo  =  9(l-9)-B0  ,  «o  =  j(l-x)-Co  .  (11) 
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Here  multipliers  at  the  arbitrary  non-negative  functions  were  taken  in  the  present  form  in 
order  to  keep  conditions  0  <  c  <  1,  0  <  0  <  1,  and  0  <  j  <  1,  when  solving  numerically  the 
constitutive  equations  of  the  model. 
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Similarly  to  [24]  we  specify  Bo  in  the  following  form: 

Bo  =  «OA-[0  -  0o(p)]  ,  (12) 

where  6o(p)  is  determined  as  an  inversed  function  of  the  equilibrium  curve  Peqv  =  Peqv(0)  and  Bo 
is  equal  to  zero  at  0  <  0o(p).  An  example  of  derivation  of  the  equilibrium  curve  has  been  shown 
in  [24]  for  sand  and  porous  2024- Aluminium.  Thus,  the  dependence  0o(p)  is  a  stationary  point 
of  the  solution  of  the  constitutive  equation  (9)  for  the  volume  concentration  parameter.  The 
magnitude  of  the  rate  ip  is  selected  in  such  a  way  that  the  pore  collapse  time  (the  relaxation 
time  for  equation  (9))  would  be  of  the  order  of  10-100  ns  as  was  assessed  in  [36].  To  satisfy  this 
assessment  for  the  relaxation  function  ip  the  constant  uoa  is  selected  as  uoa  =  2.5 xlO4  psec4.  Ao 
and  Co  are  taken  to  be  constant  because  data  for  the  mass  transformation  and  heat  exchange 
are  not  available.  As  in  the  case  of  porous  aluminium  we  select  Ao  =  0.  However,  keeping  in 
mind  the  extreme  pressure  range,  we  allow  for  a  slow  heat  exchange  with  Co  =  10-5  psecA  in 
this  case  the  rate  co  and  corresponding  rate  for  the  heat  exchange  is  of  the  order  of  1  msec  as 
assessed  in  [37], 


0  0.6  1.2  1.8  2.4  P, 

GPa 


Fig.  3.  Analytical  fit  of  the  equilibrium  curve  to  the  experimental  points  [15] 

The  data  on  the  equilibrium  curves  can  be  obtained  only  from  static  compression  data  or 
assessed  from  the  impact  data  recording  complete  stress  profiles.  Unfortunately,  these  data 
are  not  readily  available  for  the  copper  powder  with  porosity  m  =  4.  Therefore,  we  draw  an 
analytical  fit  to  the  experimental  data  [15]  for  a  porous  copper  with  m  =  1.4;  these  experiments 
are  shown  in  Fig.  3  by  solid  circles.  The  porosity  m  is  related  to  the  following  volume 
concentration  of  the  gaseous  phase  0  =  1  -  1/m  (in  this  formula  only  we  neglect  mass  of  gas 
when  comparing  it  with  mass  of  the  solid  phase  within  a  volume).  Therefore,  the 
experimental  points  in  [15]  were  found  within  the  range  of  0  from  0  up  to  1  -  (1.4)  _1  and  an 
analytical  formula  for  the  equilibrium  curve  can  be  built  up  within  this  range.  In  order  to  use 
it  for  our  case  of  the  highly  porous  copper  we  extrapolate  the  analytical  formula  up  to 
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6  =  1-  1/4=0 .75.  As  a  result,  we  pre-determine  the  analytical  function  to  be  constructed  as 
follows:  9o(p)  =  0.01215  at  p  >  2.5  GPa,  6o(p)  =  0.75  at  p  <  0.05  GPa,  and  otherwise 

6o{p)  =  avp 2  +  ci2f  +  fl3  +  ai/p  , 

where  a-\  =  -0.0269585  GPa-2,  ai  =  0.14122  GPa-1,  U3  =  0.763143,  and  «4  =  0.161493  GPa.  The  solid 
curve  in  Fig.  1  shows  the  constructed  analytical  function  6o(p)  along  with  its  comparison  with 
the  experimental  data  [15]. 

In  order  to  close  the  model  we  must  select  EOSs  for  the  two  phases.  The  EOS  for  the  gaseous 
phase  was  chosen  in  the  following  ideal  gas  form 

e(0  (A  s)  = - v  expf  s/ 

P()aW  -  U  V/C- 


r  _P_ ^ 
\Poa  J 


The  EOS  for  the  solid  phase  was  taken  in  the  following  form  [24],  which  is  close  to  that  of  Mie- 
Gruneisen  [38]: 


e(2\p,s) 


f  \ 

a0 

z 

(  \ 

P 

/  X 

p 

-1 

-  — +  c  T 
T  G ’s1  o 

1  P 

exp  s/c  r1 

yPos  ) 

P 

l  Pos  ) 

V/  CVS  J 

(13) 


Material  constants  for  the  solid  (copper)  phase  were  taken  to  be  Co  =  3.94  km/ s, 
pos  =  8.93  g/ cm3,  Cvs  =  0.00045  }/ (g.grad),  a0  =  0.96,  (3  =  2.03,  p0  =  O.IMPa  (1  bar),  T0  =  293°K.  The 
standard  values  of  constants  were  chosen  for  air  at  normal  conditions  with  the  exponent 
y  =  1.4,  thermal  capacity  cva  =  0.001  J/(g-grad),  initial  density  poa  =  0.0012  g/cm3,  and  the 
constant  So  is  calculated  from  the  condition  of  normal  pressure  at  the  initial  density  of  0.1  MPa 
(1  atm);  ep  adjusts  specific  energies  of  the  phases  at  the  normal  conditions. 

The  pressure  response  of  an  isotropic  material  to  uniform  omni-lateral  dynamic  compression 
(the  (p,  V)-curves)  is  always  associated  with  a  compression  rate  defined  as  follows 


du,  du-,  du , 

£  = - L  + - ^  + - -  . 

Sxj  dx2  dx2 


(14) 


For  the  quasi-static  response  the  rate  (14)  is  close  to  zero;  in  contrast,  for  the  shock 
compression  response  it  exceeds  thousands  of  inverse  seconds.  In  order  to  compare  the  (p, In¬ 
curves  with  experimental  pressure  response  data  available  in  the  literature,  these  curves  have 
been  calculated  with  the  system  of  equations  reduced  to  the  case  of  uniform  omni-lateral 
compression.  Equations  for  this  loading  case  are  obtained  from  the  general  system  of 
equations  with  the  assumptions  of  a  constant  strain  rate  and  uniform  spatial  deformation.  As 
a  result,  the  system  (4),  (7-9)  is  reduced  to  the  following: 
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d  In  p  _  .  dc 

^  ? 

dt  dt 


~<P* 


dd_ 

dt 


-V , 


^L  =  -co  —=  Cc(p  +  e(,y/  +  6xCO 
dt  dt  T 


(15) 


0.1  0.12  0.14  V,  cm3/g 


1  1.2  m 


Fig.  4.  The  compression  curves  and  comparison  with  the  experimental  data  [14-16]  (solid  points) 

Results  of  calculations  using  the  system  (15)  at  the  strain  rate  of  103  s-1  are  shown  in  Fig.  4, 
where  specific  volume  is  defined  as  V  =  1/p.  As  it  is  seen  from  (15),  solution  of  the  system  at  a 
constant  £  describes  the  evolution  of  parameters  with  increments  of  specific  volume 
following  increments  of  time.  Initial  data  for  parameters  c  and  0  are  determined  from  initial 
porosity,  y  is  taken  initially  to  be  0.5  in  order  to  equally  allow  for  the  heat  exchange  between 
the  phases,  and  s  =  0  (corresponds  to  normal  initial  temperature).  The  functions  and  constants 
Ao,  Bo,  and  Co  were  selected  as  above  with  the  parameter  aoa  in  (12)  to  equal  2.5><104psec-1. 
Calculated  (p,V)- curves  in  Fig.  4(a)  correspond  to  initial  densities  poo  =  6.82  g/cm3  (curve  1) 
and  poo  =  6.052  g/cm3  (curve  2)  for  porous  copper.  The  curves  are  compared  with 
experimental  data  (the  solid  points).  The  solid  circles  in  Fig.  4(a)  are  data  from  [14]  for  a 
porous  copper  with  initial  density  poo  =  6.052  g / cm3;  it  should  be  noted  that  those  data  are 
related  to  pressure  behind  the  shock  front  of  the  main  wave  versus  'final'  density  measured 
from  the  high  velocity  impact  experiments  [14-15],  using  a  gas  gun.  The  square  points  are  data 
from  [16]  obtained  in  a  similar  way  for  a  porous  copper  with  initial  densities  between  6.31 
and  6.82  g / cm3.  It  should  be  noted  that  the  loading  conditions  for  these  data  can  only  be  quite 
roughly  approximated  by  a  uniform  omni-lateral  load  with  a  constant  strain  rate. 
Nevertheless,  the  compression  curves  are  quite  close  to  the  experimental  data  seen  in  Fig.  4. 


11 


DSTO-TR-2026 


For  relatively  low  velocities  of  impact,  we  believe  that  this  is  a  reasonable  fit  for  the  case  of  a 
dynamic  low  strain  rate  compression. 

The  calculated  curve  in  Fig.  4(b)  shows  pressure  against  current  porosity  for  the  porous 
(under-compacted)  copper  with  initial  porosity  m  =  1.38  (poo  =  6.43  g/cm3).  Experimental 
points  in  Fig.  4(b)  for  a  porous  copper  with  initial  density  of  6.43  g/ cm3  have  also  been 
obtained  from  gas  gun  experiments  [15]  and  the  same  points  were  used  in  order  to  fit  the 
curve  shown  in  Fig.  3.  This  kinetic,  which  was  obtained  by  fitting  the  constitutive  equations 
within  a  range  of  porosities  close  to  1.4,  will  also  be  used  for  modelling  the  behaviour  of  shock 
compression  of  a  highly  compressed  porous  copper  with  porosity  m  =  4. 

It  should  be  noted  that  the  curves  in  Fig.  4(a)  quickly  converge  at  higher  pressures  to  the 
compression  curve  corresponding  to  the  solid  copper  regardless  of  initial  porosity  and  strain 
rate.  It  is  interesting  that  the  (p,  V)- curves  are  non-convex,  however,  the  equilibrium  adiabatic 
condition  d2p/8V2  >  0  at  s  =  const  takes  place  [24]  and  no  singularities  of  the  equation  of  state 
are  required  to  describe  this  response.  It  should  be  noted  that  a  phase  transition  could  be 
described  in  the  same  way  as  well  (see,  for  instance,  a  model  considered  in  [39]). 


5.  Modelling  of  the  shock  behaviour  of  porous  copper 

The  present  numerical  analysis  is  conducted  within  a  one-dimensional  set-up,  using  the 
Godunov  method  [40].  EOS  and  constitutive  relations  are  used  in  the  form  specified  in  the 
previous  Section.  The  Riemann  problem  (a  necessary  step  of  the  Godunov  method)  for  a 
porous  material  is  solved  by  the  method  published  in  [41],  which  is  briefly  described  in 
Appendix  A. 

First,  we  analyse  the  wave  structure  in  a  porous  copper.  For  a  pre-defined  equilibrium  curve 
the  compaction  time  in  the  shock  wave  is  affected  by  the  time  of  pore  collapse  that  is 
associated  with  deformation  of  the  material  matrix.  In  turn,  the  compaction  time  is  affected  by 
the  parameter  aoa.  Therefore,  we  calculate  the  wave  propagation  in  porous  copper  with 
parameters  and  kinetics  selected  in  the  previous  Section.  The  only  variable  will  be  the 
parameter  aoa.  The  initial  density  corresponds  to  the  porosity  of  m  =  4  (poo  =  2.232  g/ cm3). 

In  order  to  observe  how  the  pore  collapse  time,  that  is  determined  by  the  compression  kinetic, 
affects  the  wave  structure,  we  draw  the  wave  profiles  in  Fig.  5  for  the  case  corresponding  to  a 
relatively  slow  compaction  («oa  =  2. 5  x  101  psec-1) .  Consecutive  velocity  and  pressure  profiles  in 
a  sample  (the  sample's  thickness  is  10  mm)  of  the  material  being  compressed  by  a  rigid  plate 
moving  with  a  constant  velocity  Uo  correspond  to  the  intervals  shown  in  psec  against  each 
profile  throughout  the  present  section. 
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0  x,  cm  10  x,  cm  1 


Fig.  5.  The  velocity  (a)  and  stress  (b)  profiles  at  Uo  =  0.3  km/sfor  the  slow  compression  kinetics 

Similar  calculation  has  been  conducted  for  a  material  with  a  quicker  response  to  compression 
(flOA  =  2.5><104psec-1),  which  corresponds  to  a  more  realistic  pore  collapse  time  as  mentioned  in 
the  previous  Section.  Results  are  shown  in  Fig.  6  and  demonstrate  a  significant  distinction  in 
the  shock  front  width;  however,  the  character  of  the  loading,  including  evolution  of  the 
precursor  and  pressure  attenuation  of  the  main  wave,  has  not  changed  significantly.  It  should 
be  noted  that  variable  x  for  the  wave  profile  diagrams  is  the  Eulerian  coordinate  stretched  to 
the  initial  1  cm  space  for  every  profile. 


0  x,  cm  10  x,  cm  1 


Fig.  6.  The  velocity  (a)  and  stress  (b)  profiles  at  Uo  =  0.3  km/sfor  the  fast  compression  kinetics 

It  should  be  noted  that  the  slow  kinetic  provides  a  pore  collapse  time  interval  of  the  order  of 
10  psec,  whereas  for  the  fast  kinetic  the  collapse  is  a  thousand  times  faster.  This  collapse  time 
is  related  to  the  material  structure;  in  particular,  to  the  particle  size  and  it  is  comparable  with 
the  closure  time  of  a  pore.  For  example,  for  a  50pm-porous  material  loaded  by  a  shock  wave 
with  the  pore  walls  moving  with  velocity  of  the  order  1  km/ s  the  closure  time  is  0.05  psec. 
This  assessment  is  in  agreement  with  the  pore  collapse  times  derived  in  [36, 37]  as  mentioned 
earlier.  Therefore,  the  fast  kinetic  seems  to  be  the  more  appropriate  one. 
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In  order  to  compare  calculated  shock  profiles  with  experimental  data  we  use  information 
obtained  in  papers  [14-16],  Unfortunately,  the  observed  gauge  traces  are  not  available  but  the 
authors  provided  a  schematic  of  the  wave  profiles  with  several  control  points  of  the  profile 
contour.  The  gas  gun  experiments  [14]  were  conducted  with  a  plane  sandwich  set-up:  a  quartz 
6.35  mm-flyer  plate  collided  with  porous  copper  (poo  =  6.43  g/ cm3)  samples  backed  by  another 
6.35  mm  quartz  plate.  The  stress  profiles  were  recorded  at  the  two  quartz  sample  interfaces. 
Results  from  the  gauges  (the  control  points  of  the  wave  profile  contours  for  tests  RA5  and  RA2 
from  [14])  are  shown  in  Fig.  7  by  solid  squares  (test  RA2  with  the  impact  velocity  of  0.1845 
km/  s  and  sample  thickness  0.83  mm)  and  by  solid  circles  (test  RA5  with  the  impact  velocity  of 
0.3073  km/ s  and  sample  thickness  1.03  mm).  The  points  located  at  approximately  0.1  psec  in 
Fig.  7  indicate  the  experimental  pressure  levels  taken  by  the  first  gauge  after  equilibrium  on 
the  quartz  sample  surface  (the  impact  surface).  The  remaining  points  relate  to  the  wave 
recorded  by  the  second  gauge  from  the  sample  quartz  interface  of  the  target  assembly. 


Fig.  7.  The  stress  profiles  of  the  test  set-ups  compared  with  the  control  points  from  the  gas  gun 
experiments  [14] 


In  order  to  conduct  calculations,  the  solid  components  of  the  assembly  should  be  simulated  as 
well.  For  the  solid  materials  (components  of  the  flyer  plate,  screen,  backing  plate,  gauge,  etc) 
the  visco-elastic  rate  sensitive  model  [42]  has  been  employed.  The  EOS  [38]  was  used  in  the 
present  work  in  the  following  form  (which  is  close  to  (13)): 


E(p,D,s) 


2a; 


•  (Sa'  -l)2+2 b20SAD+cvT0Sn 


1 


(16) 
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Here  D  is  a  second  invariant  of  the  elastic  strain  tensor  [38]  and  5  =  p/po,  where  po  is  initial 
density.  The  strength  effects  are  described  by  the  model  with  use  of  the  relaxation  time 
function  of  shear  stresses.  A  procedure  of  how  this  function  can  be  obtained  using  two 
different  yield  limit  points  at  two  different  strain  rates  has  been  described  in  [43].  Thus,  the 
following  material  constants  are  required  by  the  model  for  description  of  the  behaviour  of 
solid  materials:  density  po,  bulk  modulus  K  =  po  Co2,  shear  modulus  G  =  po  bo2,  the  thermal 
capacity  coefficient  cv,  two  compressibility  coefficients  (exponents  in  (16))  ai  and  (3i  for  the 
modules  K  and  G,  and  the  Gruneisen  coefficient  yi.  In  the  case  of  quartz  components  we 
selected  constants  in  the  relaxation  function  in  such  a  way  that  the  relaxation  time  would  be 
large  enough  during  the  deformation  process;  then,  the  quartz  responds  nearly  elastically 
when  interacting  with  the  porous  samples. 

For  quartz  we  selected  the  following  constants  for  the  EOS  (16);  po  =  2.203  g/ cm3,  K  =  37  GPa, 
G  =  32  GPa,  cv  =  0.01  J/  (g-grad),  cu  =  1.07,  [3i  =  3,  and  yi  =  1.117.  We  used  the  fast  compaction 
kinetic  («oa  =  2.5 xlO4  psec-1)  in  all  calculations  below.  The  wave  structure  tracked  by  the 
second  gauge  at  the  rear  side  of  the  sample  gives  the  start  of  the  precursor  (solid  points  nearly 
0.3  psec  in  Fig.  7),  pressure  level  of  the  second  step  behind  the  precursor  (time  of  arrival  and 
pressure  level  identified  by  the  solid  points  at  the  breaks  of  dashed  lines),  the  middle-point  of 
the  main  compression  wave  (solid  points  at  nearly  1.3  psec  in  Fig.  7),  and  amplitude  of  the 
main  wave  (indicated  by  solid  points  at  t  =  1.8  psec).  Thus,  the  experimental  profiles  are 
characterised  schematically  by  the  dashed  lines  passing  through  the  solid  circles  for  the  test 
RA2  and  through  the  solid  squares  for  the  test  RA5. 

Corresponding  calculations  are  drawn  by  line  1  for  the  test  RA2  and  by  line  2  for  the  test  RA5. 
It  is  remarkable  that  the  second  pressure  increase  after  the  precursor  is  observed  in  the 
calculation  as  well;  it  appears  from  the  calculations  that  this  second  rise  is  actually  associated 
with  the  pressure  build  up  due  to  interaction  of  the  reflected  precursor  from  the  sample- 
quartz  interface  with  the  main  compression  wave  and  transmission  of  the  second  double-load 
precursor  to  the  interface  after  the  interaction.  It  is  interesting  that  similar  precursor 
reverberation  has  been  observed  in  reverberation  tests  [16].  Distinct  precursors  have  been 
observed  experimentally  in  other  porous  materials  as  well  (e.g.,  for  aluminium  powder, 
typical  wave  profiles  have  been  reported  in  [23]).  Thus,  we  believe  that  the  constants  selected 
for  porous  copper  are  satisfactory  for  description  of  these  features  of  the  powder  shock 
compression  and  all  subsequent  analysis  will  be  conducted  with  the  present  constants. 

Further  analysing  profiles  obtained  in  Fig.  6,  we  can  draw  the  density  and  volume 
concentration  profiles  as  shown  in  Fig.  8.  It  is  seen  that  for  the  case  of  a  relatively  low  impact 
load  the  wave  non-equilibrium  is  chiefly  governed  by  interaction  of  the  precursor  wave 
reflected  from  the  free  surface  with  the  main  compression  wave.  When  evolving,  the  particle 
velocity  behind  the  front  gradually  rises  (after  interaction  with  the  precursor),  pressure  drops, 
and  density  remains  unchanged. 
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Fig.  8.  The  density  (a)  and  volume  concentration  (b)  profiles  for  the  wave  at  Uo  =  0.3  km/s 


When  increasing  the  load  strength  the  wave  structure  changes.  The  following  two  figures 
show  profiles  for  the  waves  when  loaded  with  a  driver  velocity  of  Uo  =  1  km/  s.  Consecutive 
wave  profiles  in  Figs.  9-10  correspond  to  1.5,  3,  4.5,  and  6  psec. 


1.1 
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Fig.  9.  The  velocity  (a)  and  stress  (b)  profiles  for  the  wave  at  Uo=  1  km/s 

The  influence  of  the  precursor  in  this  case  is  still  noticeable  but  not  so  strong  because  the 
precursor  velocity  is  only  double  the  velocity  of  the  main  compression  wave.  In  this  case 
evolution  of  the  wave  profile  is  almost  entirely  governed  by  the  compaction  process  and  this 
only  slightly  alters  the  pressure  profile  from  uniformity.  Corresponding  density  and  volume 
concentration  profiles  are  drawn  in  Fig.  10.  It  is  seen  from  Fig.  9  that  the  particle  velocity  is  not 
significantly  changing  within  any  particle  in  this  nearly  Lagrangian  coordinate  system  (the 
apparent  non-uniformity  in  the  space  between  the  wave  profiles  should  be  commented  as  this 
is  associated  with  the  coordinate  transformation;  in  fact,  after  reduction  of  the  profiles  to  the 
real  Eulerian  position,  the  distance  between  profiles  is  practically  constant  with  time). 
However,  pressure  gradually  drops,  allowing  for  rearrangement  of  the  solid  matrix  and 
equilibrating  pressure  within  the  matrix  and  air.  At  the  same  time,  density  is  gradually 
increasing  along  with  the  volume  concentration  of  the  solid  phase  due  to  expansion  of  the 
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initially  compressed  solid  phase  into  the  gaseous  phase.  From  the  viewpoint  of  state  evolution 
behind  the  shock  front,  the  velocity  and  pressure  drop,  whereas  density  remains  unchanged. 
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Fig.  10.  The  density  (a)  and  volume  concentration  (b)  profiles  for  the  wave  at  Uo  =  1  km/s 

Pressure  is  equilibrating  relatively  quickly  via  relaxation  of  the  parameter  0  where  a 
stationary  point  is  achieved  at  ee  =  0  (i.e.,  when  the  phase  pressures  equilibrate).  Within  this 
range  of  loads  the  material  state  behind  the  front  undergoes  a  slow  heat  exchange, 
equilibrating  the  temperature  via  relaxation  of  the  parameter  y  with  a  stationary  point 
achieved  at  ex  =  0  (i.e.,  when  the  phase  temperatures  equilibrate).  The  parameter  evolution 
directly  behind  the  main  shock  front  might  be  explained  by  the  two-stage  character  of  the  load 
headed  by  the  precursor.  Comparing  two  elements  of  the  material,  the  first  of  which  is  close  to 
the  loading  surface  and  the  second  is  in  the  middle  of  the  sample  we  can  analyse  the  process 
of  deformation.  Both  elements  are  compressed  to  approximately  the  same  volume.  However, 
the  first  is  loaded  by  a  (practically)  one-step  load;  therefore,  the  pressure  and  temperature  are 
higher  then  for  the  second  element,  for  which  the  loading  is  closer  to  isentropic  (farther  from 
the  shock-like).  The  second  element  is  subjected  to  a  two-step  loading  because  the  relaxation  is 
taking  a  longer  time.  These  processes  are  manifested  in  the  whole  response  of  material  shown 
in  the  calculations  of  Figs.  9  and  10.  It  should  be  noted  that  the  deformation  physics  is 
associated  very  closely  with  the  mechanisms  of  compression  and  heat  exchange  governed  by 
the  kinetic  chosen.  Another  example  of  the  wave  behaviour  will  be  given  in  the  Discussion 
Section. 

To  complete  the  analysis,  propagation  of  a  fairly  strong  shock  wave  has  also  been  calculated 
with  a  driving  plate  velocity  of  Uo  =  3  km/ s.  Consecutive  profiles  in  Figs.  11-12  correspond  to 
the  time  moments  of  0.5,  1,  1.5,  and  2  psec.  It  is  seen  from  Fig.  11  that  in  this  case  the 
precursor  effect  is  almost  negligible.  Corresponding  profiles  for  density  and  volume 
concentration  of  the  gaseous  phase  are  drawn  in  Fig.  12.  It  is  also  seen  from  these  figures  that 
in  this  case  the  shock  wave  is  almost  fully  stationary  with  a  relatively  small  drift  of  velocity 
and  pressure  directly  behind  the  shock  front  due  to  the  heat  exchange  effects.  It  should  be 
noted  though  that  the  trend  remains  the  same  as  in  the  preceding  calculation.  If  the  load  is 
further  increased,  the  deformation  is  closer  to  a  one-step  loading  because  the  main  wave 
overtakes  the  precursor. 
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Fig.  11.  The  velocity  (a)  and  stress  (b)  profiles  for  the  wave  at  UO  =  3  km/s 


Thus,  the  constitutive  modelling  of  shock  propagation  in  porous  copper  shows  a 
disagreement  between  the  parameters  (Ul,  Pl,  Pl)  at  the  loading  surface  (piston)  and  the 
parameters  (Uf,  Pf,  Pf)  directly  behind  the  shock  front.  The  characteristic  points  of  the  profile  4 
in  Figs.  11-12  are  marked  by  arrows  and  labelled  by  letters  'L'  and  'F',  correspondingly. 
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Fig.  12.  The  density  (a)  and  volume  concentration  (b)  profiles  for  the  wave  at  Uo  =  3  km/s 

It  is  interesting  to  note  that  the  percentage  deviation  (Ul  -  Uf)/Ul  increases  from 
approximately  -7%  first  (it  is  negative  in  this  case  because  of  the  precursor  effect)  to  5-7%  as 
seen  from  comparison  of  the  results  in  Figs.  6-9(a);  after  that  it  decreases  down  to 
approximately  2-3%  when  comparing  the  results  from  Fig.  9(a)  and  Fig.  11(a).  This  situation 
contrasts  with  that  for  a  solid  material  where  the  jump-like  approximation  of  the  shock  wave 
due  to  a  constant  load  is  very  good  for  the  load  range  where  the  strength  and  phase  transition 
effects  can  be  ignored.  This  is  an  alarming  fact  because  traditionally,  when  using  the  jump-like 
approximation,  the  Hugoniot  density  is  calculated  from  the  following  relation  from  system 
(1):  p  =  po  D/(D  -  U),  where  U  is  presumed  to  be  a  constant  velocity  (the  loading  velocity) 
behind  the  shock  wave.  Because  this  is  not  so  in  the  present  case,  we  evaluate  Hugoniots  in 
the  next  Section  from  direct  calculations  of  states  behind  the  shock  wave,  using  the  set-ups 
and  input  data  recovered  as  closely  as  possible  to  the  ones  published  in  the  literature. 
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6.  Brief  description  of  available  experiments 

Focusing  on  the  tests  for  porous  copper  with  porosity  m  =  4,  we  list  all  available  experimental 
data  in  Table  1.  We  limit  our  consideration  to  data  for  the  pressure  range  lower  than  200  GPa. 
The  Table  lists  17  data  points  from  publications  [6-9,  11-12],  The  Table  includes  the  IM- 
reduced  data  of  the  particle  velocity  U  and  the  shock  velocity  D  (in  km/s)  for  copper 
powders.  As  was  mentioned  earlier,  the  data  from  [6-9]  were  reportedly  subject  to  another 
reduction  that  included  a  data  correction  via  elimination  of  the  wave  attenuation  effect.  As 
was  described  in  Section  2,  the  particle  velocity  U  was  calculated  for  the  reduced  data,  using 
the  Impedance  Mismatch  (IM)  method.  The  procedure  of  calculation  of  the  particle  velocity 
with  the  IM-method  is  described  in  Appendix  B;  formula  (24)  in  the  Appendix  gives  the 
method  of  calculation  for  the  particle  velocity  Uh  =/(D,  Uo)  that  is  denoted  by  U  in  the  present 
Section.  While  using  an  inverse  procedure,  we  can  recover  the  data  Uo  in  the  screen  (from  the 
values  of  Uh  and  D  that  are  listed  in  the  experimental  works),  using  tabulated  Hugoniots  for 
the  screen  materials,  where  the  Hugoniots  are  usually  presented  as  linear  particle  velocity- 
shock  velocity  dependencies 


D  =  a  +  b  U  .  (17) 

The  screen  material  that  was  used  for  the  tests  in  [6,  7]  is  the  Russian  analogue  (AD-1)  of  the 
aluminium  1100-0  alloy  as  reported  in  [44].  The  Hugoniot  that  was  used  for  the  IM-reduction 
was  reported  in  [7]  as  follows: 

D  =  5.333+1.356  U  .  (18) 

Thus,  using  the  Hugoniot  for  drawing  the  shock  adiabat  in  the  (P,U)-space  with  the  help  of 
the  diagram  in  Fig.  2(b),  we  can  recover  state  '[2]'  that  is  characterised  by  the  parameters 
(Po,Uo).  We  will  denote  the  particle  velocity  behind  the  shock  front  by  Uscr  in  the  present 
Section.  Then,  using  USCr,  we  can  calculate  Dscr  from  (17). 

After  doing  this  with  the  assumption  that  the  screen  for  all  the  tests  is  made  from  the 
aluminium  alloy,  all  recovered  screen  data  have  been  listed  in  Table  1.  Reported  screen  data 
(the  particle  velocity)  is  available  for  point  5  only,  and  the  data  restored  from  the  (U,  D)-data 
for  this  point  reproduce  the  experimental  USCr  value  quite  well.  Screen  materials  are  not 
known  for  points  7,  11,  13,  and  17  and,  therefore,  we  assumed  the  aluminium  alloy  for  the 
screen  by  default.  Initial  densities  for  porous  powders  of  points  7  and  11  slightly  deviate  from 
the  porosity  m  =  4;  when  restoring  the  screen  data  and  conducting  constitutive  modelling  we 
always  used  the  actual  density  values  for  the  corresponding  experiments  (these  values  are 
shown  in  the  references/ comments  column  followed  by  'p').  Points  14  and  15  were  actually 
obtained  with  iron  screens  and  the  particle  velocities  that  were  measured  in  the  experiments 
along  with  the  shock  velocities  in  the  screens  are  shown  in  the  references/ comments  column 
followed  by  'Fe'.  The  listed  screen  data  USCr  and  DSCrfor  these  points  is  a  result  of  recalculation 
for  the  aluminium  screen,  using  the  (U,  D)-pairs  reported  in  [8,  9]. 
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Table  1.  The  Hugoniot  velocity  data  from  available  experiments 


N 

User 

U 

D 

Dscr 

P 

P 

Ref.,  Comments 

1 

0.37765 

0.661 

0.946 

5.845 

1.395 

7.402 

[6] 

2 

0.765 

1.235 

1.66 

6.37 

4.572 

8.71 

[61 

3 

1.135 

1.725 

2.33 

6.87 

8.963 

8.588 

[6] 

4 

1.46 

2.12 

2.93 

7.31 

13.852 

8.066 

[6] 

5 

1.95 

2.68 

3.84 

7.97 

22.95 

7.382 

[7] 

6 

2.12 

2.89 

4.07 

8.209 

26.23 

7.6916 

[61 

7 

2.514 

3.33 

4.662 

8.74 

35.132 

7.921 

[11],  p(2.263) 

8 

2.69 

3.51 

5.10 

8.98 

39.92 

7.153 

[61 

9 

2.918 

3.77 

5.41 

9.29 

45.483 

7.356 

[61 

10 

3.288 

4.16 

6.07 

9.79 

56.31 

7.087 

[6] 

11 

3.472 

4.4 

6.248 

10.04 

60.536 

7.445 

[11],  p(2.202) 

12 

3.716 

4.61 

6.81 

10.37 

70.01 

6.903 

[6] 

13 

4.94 

6.00 

8.325 

12.033 

111.4 

8.0 

[12] 

14 

5.147 

6.15 

8.96 

12.31 

122.88 

7.1106 

[8],  Fe(4.15,11.05) 

15 

5.315 

6.42 

8.79 

12.54 

125.8 

8.27 

[9],  Fe(4.26,10.67) 

16 

5.615 

6.68 

9.51 

12.946 

141.665 

7.4937 

[6] 

17 

6.8 

8.0 

11.028 

14.56 

196.7 

8.146 

[12] 

Using  the  velocities  U  and  D  (in  km/s),  pressure  (in  GPa)  and  density  (in  g/cm3)  are 
calculated  with  the  relations  (1)  and  are  also  shown  in  Table  1  in  the  P-  and  p-columns.  Points 
1  and  2  were  possibly  measured  by  piezoceramic  sensors,  according  to  [6],  and  the  others  by 
the  pin  technique.  According  to  [6],  one  of  the  points  was  possibly  found  for  a  sample  with 
evacuated  air.  This  sample  had  to  be  cased  in  order  to  allow  for  the  evacuation.  However,  no 
information  about  the  test  conditions,  including  set-up,  initial  parameters,  and  quality 
(homogeneity)  of  the  material  after  the  air  evacuation,  was  reported,  therefore,  this  point 
cannot  be  adequately  modelled. 

It  is  seen  from  the  (P,  p)-data  in  Table  1  that  the  apparent  abnormality  is  clearly  exhibited  by 
elevated  densities  for  points  2, 3,  and  4;  all  these  points  were  reported  in  [6].  At  the  same  time, 
as  we  observed  in  the  previous  Section,  the  velocity  range  covering  these  points  is  within  the 
range  where  the  velocity  non-equilibrium  behind  the  shock  front  is  most  noticeable.  This  fact 
will  be  paid  special  attention  in  the  next  modelling  Section. 

The  experimental  design  for  tests  [6]  was  described  in  greater  detail  in  recent  publications 
[1, 25].  These  papers  show  a  general  set-up  for  the  plate  impact  experiment,  which  is  close  to 
the  schematic  in  Fig.  1;  in  reality,  the  set-up  is  complemented  with  a  Plane  Wave  Generator 
(PWG)  spaced  by  an  air  gap  to  the  flyer  plate.  The  PWG  and  air  gap  dimensions  were  not 
reported.  The  screen  material  for  the  majority  of  experiments  was  reported  to  be  the 
aluminium  1100-0  alloy;  two  points  14  and  15  were  reported  in  [8, 9]  to  be  obtained  with  steel 
screens.  Because  the  abnormality  was  observed  for  tests  with  the  aluminium  screens,  we  pay 
special  attention  to  the  set-up  with  the  aluminium  screen.  According  to  [6, 25-27],  the  sample 
thickness  varied  from  1.5  to  5  mm  and  for  the  flyer  and  screen  plates  varied  from  2  to  5  mm. 
Keeping  in  mind  that  the  shock  attenuation  effects  are  minimal  for  the  thickest  flyer  and 
screen  plates  and  the  thickness  distribution  over  the  experiments'  set  is  not  reported,  we 
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selected  a  set-up  with  favourable  geometrical  parameters  in  order  to  minimise  the  wave 
attenuation,  i.e.,  thickness  of  the  plates  was  selected  to  be  5  mm.  It  is  expected  that  thickest 
sample  correspond  closer  to  stationary  conditions,  therefore,  for  the  basic  set-up  we  selected 
the  sample  thickness  to  be  5  mm  as  well. 

In  order  to  analyse  possible  effects  of  the  set-up  geometry  variation  we  set  the  screen 
thickness  to  be  5  mm  and  selected  two  options  for  the  flyer  plate  thickness.  The  first,  where 
the  wave  attenuation  effects  are  possible,  is  with  a  flyer  plate  thickness  of  5  mm.  The  second  is 
with  thicker  flyer  plates  (from  6  to  15  mm,  depending  on  the  impact  velocity),  which 
eliminates  the  wave  attenuation  factor. 


7.  Results  of  modelling  the  shock  compression  tests 

This  Section  simulates  the  conditions  of  the  shock  compression  tests  described  in  the  previous 
Section  as  closely  as  possible.  Because  the  specifications  of  the  explosive  charges  accelerating 
the  flyer  plates  were  not  known,  we  assumed  that  each  flyer  plate  had  a  constant  velocity 
2Uscr  with  User  specified  in  respective  column  of  Table  1.  Two  cases  were  considered  for  the 
study  of  possible  influences  of  the  wave  attenuation  on  the  calculation  results.  Conditions  of 
the  first  case  presume  that  thickness  of  the  flyer  plate  is  variable  such  that  the  time  h,  when 
the  shock  wave  in  the  porous  sample  is  overtaken  by  the  rarefaction  wave  from  the  free 
surface  of  the  flyer  plate,  would  be  not  less  than  the  time  of  registration  of  a  signal  at  t  =  h  (see 
Fig.  2(a)).  For  example,  with  the  basic  set-up  for  low  impact  velocities  (points  1-5  in  Table  1) 
the  flyer  plate  thickness  should  be  not  less  than  15  mm  and  for  high  impact  velocities  (points 
13-15)  -  not  less  than  6  mm.  The  second  case  reproduces  the  best-reported  experimental  set-up 
(from  the  viewpoint  of  minimal  wave  attenuation)  with  a  fixed  flyer  plate  thickness  of  Wf  =  5 
mm  (Fig.  2(a)).  The  screen's  thickness  was  fixed  to  Wa  =  5  mm.  For  the  basic  set-up  the 
sample's  thickness  was  selected  to  be  5  mm  as  well  (Ws  =  5  mm  in  Fig.  2(a)). 

Determination  of  the  shock  wave  arrival  numerically  is  subjective  because  it  is  associated  with 
a  preferred  choice  of  the  characteristic  point  on  the  shock  profile  in  order  to  calculate  the 
profile's  velocity.  To  eliminate  this  factor,  we  attempted  to  reproduce  the  experimental  set-up 
down  to  the  level  of  the  gauge  assembly.  To  realise  this  in  the  calculation,  we  attached  to  the 
rear  surface  of  a  sample  (at  x  =  Wf  +  Wa  +  Ws  in  Fig.  2(a))  a  thin  aluminium  plate  of  thickness 
Wb  =  0.2  mm  (see  Fig.  2(a)),  which  could  be  associated  with  the  end  wall  of  the  protective  cup 
E  of  the  gauge  G  in  Fig.  1.  In  the  calculation  runs  we  traced  the  free  rear  surface  of  the  plate 
until  it  reached  the  point  x  =  Wf  +  Wa  +  Ws  +  Wb  +  Wg  at  t  =  te,  where  the  gap  space  Wg  was 
selected  to  be  0.2  mm  as  well.  This  time  is  considered  to  be  the  time  of  arrival  of  the  shock 
wave  to  the  free  surface  of  the  sample  (at  x  =  Wf  +  Wa  +  Ws  in  Fig.  2(a)). 

Normally  two  runs  were  calculated  where  the  first  one  was  with  the  sample  removed  from 
the  numerical  set-up.  For  this  run  the  time  of  arrival  U  of  the  wave  at  the  point  x  =  Wf  +  Wa 
was  incremented  by  the  time  of  propagation  of  the  shock  wave  through  the  gauge  cup  and  the 
time  of  closing  the  gap  Wg  (for  the  set-up  with  the  sample  shown  in  Fig.  2(a)  this  increment 
corresponds  to  f6  -  h),  resulting  in  the  time  moment  fi'.  Thus,  the  averaged  shock  wave 
velocity  Dcai  in  a  porous  sample  was  calculated  as  Ws/(  te  -  h  ')■  The  material  states  for 
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calculation  of  Hugoniots  were  taken  at  the  time  moment  when  the  shock  wave  starts  moving 
the  sample-cup  interface  (an  element  of  the  material  at  x  =  Wf  +  Wa  +  Ws).  Similarly,  the  shock 
velocity  in  the  screen  was  calculated  in  the  same  set-up  by  simply  replacing  the  investigated 
material  T  in  the  area  Wf  +  Wa  <  x  <  Wf  +  Wa  +  Ws  by  the  standard  material  S.  It  should  be 
noted  that  in  our  calculations  we  tested  aluminium,  copper,  or  steel  as  the  cup's  material,  and 
we  did  not  notice  an  essential  difference  in  calculated  shock  velocity,  as  the  cup's  specific 
mass  did  not  change  dramatically. 

For  the  aluminium  1100-0  alloy  standard  material,  we  selected  the  following  constants  for  the 
EOS  (16);  p0  =  2.712  g/cnP,  K  =  68.8  GPa,  G  =  27  GPa,  cv  =  0.0009  J/(g-grad),  cn  =  1.0, 
Pi  =  3.577,  and  yi  =  1.088.  Two  points  used  for  the  shear  stress  relaxation  kinetics  fitted  with 
the  use  of  the  method  [44]  were  taken  from  [45]  as  the  yield  limit  Y  against  strain  rate: 
Yi  =  0.0965  GPa  at  E\  =  1.8  sec-1  and  Y2  =  0.124  GPa  at  £2  =  2600  sec-1.  The  constitutive 
equations  for  the  porous  copper  were  used  with  the  fast  compaction  kinetics,  specifying 
rtoA  =  2.5><104  psec-1. 

Table  2.  Calculation  results  for  the  basic  set-up  and  comparison  with  the  experiments 


N 

Deal 

Dexp 

Ucal 

U'c 

Uexp 

Dscai 

Dsr 

p  cal 

pred 

1 

0.995 

0.946 

0.64 

0.66 

0.661 

5.78 

5.845 

5.86 

7.402 

2 

1.837 

1.66 

1.173 

1.21 

1.235 

6.27 

6.37 

6.06 

8.71 

3 

2.55 

2.33 

1.64 

1.69 

1.725 

6.79 

6.87 

6.19 

8.588 

4 

3.14 

2.93 

2.04 

2.09 

2.12 

7.25 

7.31 

6.29 

8.066 

5 

3.96 

3.84 

2.6 

2.66 

2.68 

7.95 

7.97 

6.46 

7.382 

6 

4.24 

4.07 

2.79 

2.86 

2.89 

8.19 

8.209 

6.51 

7.6916 

7 

4.88 

4.662 

3.33 

3.3 

3.33 

8.76 

8.74 

6.68 

7.921 

8 

5.14 

5.10 

3.43 

3.5 

3.51 

9.01 

8.98 

6.69 

7.153 

9 

5.49 

5.41 

3.7 

3.75 

3.77 

9.34 

9.29 

6.76 

7.356 

10 

6.06 

6.07 

4.09 

4.16 

4.16 

9.87 

9.79 

6.87 

7.087 

11 

6.35 

6.248 

4.3 

4.39 

4.4 

10.13 

10.04 

6.871 

7.445 

12 

6.7 

6.81 

4.55 

4.63 

4.61 

10.48 

10.37 

6.99 

6.903 

13 

8.48 

8.325 

5.9 

5.97 

6.00 

12.23 

12.033 

7.32 

8.0 

14 

8.8 

8.96 

6.1 

6.18 

6.15 

12.52 

12.31 

7.37 

7.1106 

15 

9.04 

8.79 

6.29 

6.37 

6.42 

12.76 

12.54 

7.41 

8.27 

16 

9.44 

9.51 

6.6 

6.69 

6.68 

13.19 

12.946 

7.48 

7.4937 

17 

11.1 

11.028 

7.9 

7.97 

8.0 

15.06 

14.56 

7.74 

8.146 

The  state  behind  the  shock  front  was  assessed  in  the  same  manner  as  the  state  F  (will  also  be 
referred  to  as  F-state  below)  in  Figs.  11-12.  Calculation  results  for  the  first  case  (when  setting 
up  thick  flyer  plates  in  order  to  eliminate  the  wave  attenuation  effect)  are  listed  in  Table  2 
along  with  experimental  shock  velocities,  IM-reduced  particle  velocity  and  density  (pred) 
transferred  from  Table  1.  In  the  Table,  Dcai  and  Dexp  are  the  calculated  and  experimental  shock 
velocities  in  porous  sample,  Ucai  and  UexP  are  the  calculated  (F-velocity)  and  experimental  (IM- 
reduced)  particle  velocities  in  the  sample,  Dscai  and  Dsr  are  the  calculated  and  experimental 
(IM-recovered  from  Dexp  and  Uexp)  shock  velocities  in  the  screen,  and  pcai  is  the  calculated 
density  (F-density)  in  the  sample.  It  is  seen  that  the  velocities  are  in  a  good  agreement. 
Moreover,  in  order  to  eliminate  the  screen  effect  because  of  a  possible  mismatch  of  a  material 
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described  with  the  use  of  the  present  EOS,  with  the  actual  aluminium  material  used  in  the 
tests,  we  can  repeat  the  procedure  of  elimination  of  the  screen  data  by  recalculation  of  the 
particle  velocity  U'c  in  the  material  investigated,  using  the  equation  (24)  in  Appendix  B. 

Then,  using  calculated  shock  velocity  Dcai  we  can  draw  the  Rayleigh  line  (shown  in  Fig.  2(b) 
with  the  shock  velocity  denoted  by  Dc)  and  compare  the  velocity  U'c  (denoted  by  U'  in  Fig. 
2(b))  with  the  experimental  values  (these  values  were  obtained  in  [6-12]  by  the  IM-reduction 
procedure  from  un-reported  experimental  shock  wave  velocities  in  the  screen  and  in  the 
porous  sample  followed  by  a  post-correction  for  the  wave  attenuation  effects).  It  is  seen  that 
after  applying  this  screen  elimination  correction,  the  agreement  of  U'c  and  Uexp  is  remarkable 
(the  largest  disagreement  for  the  whole  dataset  is  approximately  2%).  However,  the  densities 
are  in  a  significant  disagreement  (for  the  first  6  points  the  minimum  disagreement  is  12.5% 
and  the  largest  is  over  30%).  Therefore,  the  procedure  of  density  determination  is  likely  to  be 
an  issue;  as  error  estimates  in  Appendix  B  show,  the  density  error  increases  with  the  porosity 
proportionally.  Similar  error  estimates  for  density  are  quite  common  and  reported  in  various 
sources  (e.g.,  see  [1,  21,  46]). 

The  velocity  data  from  Table  2  are  summarised  in  Fig.  13,  where  crosses  correspond  to  the 
experimental  points  Dexp  (shown  in  the  figures  below  as  Dp  for  the  tested  material  and  Ds  for 
the  screen  material)  and  Uexp;  the  solid  circles  correspond  to  the  present  calculation  results  Dcai 
and  Ucai  for  the  tested  material  and  Dscai  and  U  (not  listed  in  Table  2)  for  the  screen  material. 

10 

8 

6 

4 

2 

o 

Fig.  13.  The  (D-D)-dependence for  the  basic  set-up  (a)  and  the  (D-U)-dependencies  for  the  standard 
(b)  and  for  the  investigated  material  (c)  with  elimination  of  the  shock  attenuation  effects.  The 
calculations  (solid  points)  are  compared  with  the  experiment  (crosses) 

For  the  case,  when  thickness  of  the  flyer  plate  is  fixed  to  5  mm,  the  wave  attenuation  effects 
are  more  significant  at  low  impact  velocities  because  the  shock  velocity  Ds  in  the  porous 
material  is  substantially  lower  than  the  sound  velocity  in  the  standard,  resulting  in  the  time  of 
start  of  the  wave  attenuation  (t3  in  Fig.  2(a))  to  be  less  than  the  time  of  arrival  of  the  shock 
wave  to  the  free  surface  (f=f4).  Corresponding  calculated  data  are  shown  in  Fig.  14.  It  is  seen 
from  Fig.  14  that  the  velocity  results  are  quite  similar  to  the  previous  case,  however,  the 
influence  of  the  wave  attenuation  on  the  (D-U)-dependence  for  the  tested  material  in  the 
region  of  maximum  influence  of  the  wave  attenuation  factor  (the  low  impact  velocity  region) 
is  more  distinct. 
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Fig.  14  The  (D-D)-dependence  for  the  basic  set-up  (a)  and  the  (D-U)-dependencies  for  the  standard 
(b)  and  for  the  investigated  material  (c)  with  the  shock  attenuation  effects  in  place.  The 
calculations  (solid  points)  are  compared  with  the  experiment  (crosses) 

It  should  be  mentioned  that  the  disagreement  between  the  calculated  and  experimental  shock 
velocity  might  appear  quite  significant  (up  to  9%  in  a  worst  scenario).  However,  when 
running  calculations  with  the  wave  attenuation  effects  there  is  a  deviation  in  the  velocity 
values;  this  deviation  exhibits  an  essential  contribution  of  the  wave  attenuation  correction 
procedure,  unfortunately,  this  procedure  was  not  documented  for  the  present  experiments. 
For  example,  the  calculated  shock  velocities  for  the  first  and  second  points  (with  wave 
attenuation)  transform  to  Dcai  =  0.843  km/  s  and  Dcai  =  1.77  km/ s;  this  gives  a  variability  for  the 
first  point  from  +5%  to  -10%  and  a  variation  from  +9%  to  +6%  for  the  second  point,  and  so  on. 
Such  variations  quickly  decrease  when  moving  from  the  non-equilibrium  low  pressure  area 
into  the  area  where  the  jump-like  approximation  is  more  appropriate;  for  example,  for  the 
data  in  Table  2,  the  disagreement  decreases  to  1-2%  for  data  points  7  and  above. 

In  order  to  more  clearly  demonstrate  the  difference  between  the  two  cases  after  data  IM- 
reduction,  we  draw  the  Hugoniots  ((P,V)-response  curves).  These  curves  are  shown  in  Fig.  15. 
The  first  set  of  points  (circles  through  which  a  curve  1  is  drawn)  corresponds  to  calculated 
pressure  and  density  behind  the  shock  front  (the  F-pressure  and  F-density).  The  second  set  of 
points  (squares  through  which  a  curve  2  is  drawn)  correspond  to  pressure  and  density 
calculated  according  to  the  mass  and  momentum  balance  across  the  shock  wave  expressed  by 
the  system  of  equations  (1);  the  shock  velocity  and  particle  velocity  values  were  taken  from 
the  present  calculations  as  Dcai  for  the  shock  velocity  and  U'c  for  the  particle  velocity.  Crosses 
in  Fig.  15  show  the  experimental  (P-V)-points  (IM-reduced  data)  listed  in  the  literature. 
Fig.  15(rt)  shows  the  Hugoniots  for  the  case,  where  the  shock  attenuation  effects  are 
eliminated;  Fig.  15(b)  shows  similar  curves  for  the  basic  set-up  with  shock  attenuation 
affecting  the  data. 

The  calculated  Hugoniots  show  quite  a  significant  deviation  of  the  states  achieved  directly 
behind  the  shock  front  (curves  1)  from  the  states  calculated  with  the  help  of  equation  (24) 
using  given  Uo  and  Dcai  (curves  2  were  obtained  with  the  IM-procedure  of  data  reduction). 
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Fig.  15.  The  (P,V)-response  curves  for  the  case  of  negligible  wave  attenuation  (a)  and  for  the  case 
when  the  shock  attenuation  affects  the  data  (b) 

It  is  interesting  to  note  that  the  pressure  deviations  are  quite  small.  Thus,  the  density  obtained 
with  the  data  reduction  is  actually  the  major  factor  contributing  to  the  disagreement. 
Remarkably,  curve  2  in  Fig.  15(a)  traces  quite  well  the  well-known  quasi-static  behaviour  of 
porous  material;  at  the  same  time,  curve  2  in  Fig.  15(b)  demonstrates  clear  signs  of  the 
compression  abnormality  that  is  typically  reported  for  highly  porous  materials.  In  contrast, 
the  Flugoniot  states  achieved  in  the  material  behind  the  shock  front  (curves  1)  are  close  to  each 
other  regardless  of  the  set-up.  All  of  the  data,  including  the  instantaneous  Ffugoniots 
(curves  1),  calculated  IM-reduced  data  Flugoniots  (curves  2)  and  reported  experimental  data 
converge  to  each  other  as  the  pressure  increases,  when  the  non-equilibrium  effects  vanish  and 
the  jump-like  feature  prevails  in  the  shock  profile. 


8.  Discussion  on  the  non-equilibrium  and  energy 

exchange  kinetics 

In  the  present  Section  we  will  try  to  analyse  the  non-equilibrium  compression  response  in 
slightly  more  detail  from  the  viewpoint  of  its  influence  on  the  Flugoniot  curves. 

The  Flugoniot  parameters  for  materials  exhibiting  constitutive  behaviour  could  be  interpreted 
in  two  ways;  non-equilibrium  and  quasi-equilibrium.  The  non-equilibrium  (instantaneous) 
Flugoniot  state  could  be  defined  in  the  traditional  way  as  that  directly  behind  the  main  shock 
front  (such  as  the  F-state  considered  above).  The  quasi-equilibrium  Flugoniot  states  could  be 
defined  as  those  along  the  shock  path  at  any  time  behind  the  shock  front  and  they  are 
associated  with  the  load  application  time.  In  the  experiments  [6],  the  observed  shock  velocities 
were  stationary  regardless  of  the  sample  thickness.  Flowever,  it  should  be  mentioned  that 
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other  parameters  were  not  controlled  in  the  experiments  because  the  shock  velocity  was  the 
only  parameter  measured.  Therefore,  if  other  parameters  varied  behind  the  shock  front  along 
the  wave  propagation  path,  they  could  not  be  deduced  from  the  experimental  time  of  arrival 
data.  According  to  our  modelling  representations,  that  does  not  contradict  the  experimental 
results,  since  the  non-equilibrium  Hugoniot  state  can  vary  depending  on  how  long  the  shock 
wave  propagates  within  the  material.  From  this  viewpoint  the  Hugoniot  state  depends  on  the 
sample  thickness  and  in  order  to  specify  it  experimentally,  the  wave  profiles  must  be 
recorded.  At  the  same  time,  quasi-equilibrium  Hugoniot  states  change  slowly  (compared  to 
the  transformation  time  in  the  shock  wave)  along  the  wave  propagation  path.  The  longer  a 
wave  propagates  the  closer  the  states  are  to  the  equilibrium  Hugoniot. 

To  illustrate  that  the  shock  velocity  data  is  not  sufficient  for  a  full  description  of  the  Hugoniot 
in  the  (P,V)-space,  calculations  below  are  conducted  for  the  basic  set-up  with  a  modified  heat 
exchange  kinetics.  The  shock  attenuation  effects  are  eliminated  in  these  calculations  by 
appropriate  selection  of  the  flyer  plate  thickness.  Starting  from  the  energy  conservation  laws 
for  the  two  phases,  where  a  specific  heat  transfer  coefficient  h  can  be  introduced,  it  is  easy  to 
associate  the  rate  function  for  the  parameter  y  with  the  heat  transfer  kinetic  proportional  to 
/7T(Ti— T2)2/ (pTiT2S2).  In  terms  of  the  kinetic  for  the  parameter  y,  we  can  choose  the  following 
non-negative  function  coo  as  the  kinetic  function: 


(D0 


where 


n  =  2.2 -10  In 


In 


n  =  2.0 


In  (T/Tc)  >1 
otherwise  . 


The  specific  heat  transfer  coefficient  is  taken  to  be  h  =  104  W/ (m3  grad),  Tc  =  104  K. 
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The  calculated  data  for  the  shock  and  particle  velocities  are  summarised  in  Fig.  16 
(corresponding  data  for  the  previous  calculation  are  shown  in  Figs.  13(a)  and  13(c)).  The 
(D-U)-dependence  for  the  screen  material  is  identical  to  that  shown  in  Fig.  13(b).  Details  of  the 
data  are  reported  in  the  table  form  in  Table  3.  Only  the  aluminium-screen  data  points  and 
points  for  porosity  m  =  4  are  reported  in  this  case. 
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Fig.  16.  The  (D-D)-dependence  for  the  basic  set-up  (a)  and  the  (D-U)-dependencefor  the  investigated 
material  (b)  with  elimination  of  the  shock  attenuation  effects.  The  calculations  (solid  points) 
are  compared  with  the  experiment  (crosses). 


Table  3.  Results  of  calculation  with  the  modified  kinetics  and  comparison  with  experiment 


N 

Deal 

Dexp 

Ucal 

Uexp 

p  cal 

pred 

1 

0.905 

0.946 

0.51 

0.661 

5.69 

7.402 

2 

1.76 

1.66 

1.19 

1.235 

6.11 

8.71 

3 

2.35 

2.33 

1.52 

1.725 

6.18 

8.588 

4 

2.90 

2.93 

1.76 

2.12 

6.03 

8.066 

5 

3.69 

3.84 

2.23 

2.68 

5.79 

7.382 

6 

3.96 

4.07 

2.49 

2.89 

6.21 

7.6916 

8 

4.80 

5.10 

3.14 

3.51 

6.59 

7.153 

9 

5.13 

5.41 

3.31 

3.77 

6.67 

7.356 

10 

5.66 

6.07 

3.70 

4.16 

6.88 

7.087 

12 

6.24 

6.81 

4.13 

4.61 

6.96 

6.903 

13 

7.91 

8.325 

5.42 

6.00 

7.51 

8.0 

16 

8.82 

9.51 

6.13 

6.68 

7.77 

7.4937 

17 

10.38 

11.028 

7.32 

8.0 

8.09 

8.146 

It  is  seen  that  description  of  experimental  data  is  even  better  for  this  case  in  the  low  velocity 
region  (up  to  5%)  when  comparing  with  the  results  of  the  previous  calculation.  However,  in 
the  high-velocity  range  agreement  of  the  (D-U)-results  with  those  reported  in  the 
experimental  works  is  actually  worse  than  in  the  previous  case  (shown  in  Fig.  13(a)). 
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Fig.  17.  The  (P ,V)-response  curves  for  the  case  of  negligible  wave  attenuation  and  modified  heat 
transfer  kinetic 

It  should  be  remembered  that  the  particle  velocity  data  is  not  obtained  directly  in  the 
experiments.  Overall,  the  calculated  and  experimental  shock  velocity  dependencies  ((D-D)- 
data)  agree  well  for  both  calculations. 

It  is  interesting  to  compare  the  Hugoniots  for  both  cases.  The  Hugoniot  for  the  previous  case 
is  shown  in  Fig.  15(a)  and  a  similar  plot  for  the  present  case  is  shown  in  Fig.  17,  keeping  the 
same  notations.  Several  comments  can  be  made  concerning  the  calculated  Hugoniot.  Firstly, 
there  is  a  significant  pressure  difference  between  the  IM-reduced  data  Hugoniot  (curve  2),  the 
calculated  Hugoniot  (curve  1),  and  corresponding  experimental  points.  This  is  obviously 
caused  by  significant  non-equilibrium  behind  the  shock  front.  It  is  interesting  to  note  that  the 
present  Hugoniots  demonstrate  the  abnormality  feature  for  both  curves  in  the  conditions 
where  the  wave  attenuation  effect  is  eliminated. 

In  order  to  observe  the  non-equilibrium  behaviour  we  draw  the  shock  wave  profiles  for  quite 
strong  shock  wave  (velocity  of  the  driving  plate  is  3  km/s).  In  this  case  for  the  previous 
calculation  we  observed  negligible  non-equilibrium  effects  (Figs.  11-12).  In  the  present  case  we 
observe  a  strong  non-equilibrium  with  a  distinct  kinetic  effect  at  the  shock  front  as  shown  in 
Fig.  18.  Two  consecutive  profiles  in  Fig.  18  correspond  to  the  time  samples  at  0.5  and  1  psec. 
Obviously,  the  modified  kinetic  uses  a  boosted  (even  for  such  extreme  pressures  and 
temperatures)  heat  transfer  coefficient,  activating  easily  the  kinetic  behaviour  behind  the 
shock  front.  Nevertheless,  this  calculation  illustrates  quite  clearly  that  the  time  of  arrival 
(shock  velocity)  data  can  be  described  well  using  different  kinetics.  Therefore,  the  shock 
profiles  recorded  in  experiments  would  be  the  very  useful  information  for  description  of  the 
shock  compression  of  materials  exhibiting  constitutive  behaviour. 
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Fig.  18.  TIte  velocity  (a),  stress  (b),  and  density  (c)  profiles  for  the  wave  at  Uo  =  3  km/s 


In  turn,  consideration  of  a  complex  thermal  excitation  mechanism  tends  to  complicate  the 
constitutive  behaviour  of  a  material  and  activates  non-equilibrium  behaviour  behind  the 
shock  front.  In  this  case  several  factors  are  influential  in  the  Hugoniot  description  in  the 
pressure-density  space:  1)  position  of  the  point,  where  data  are  taken  from,  from  the  loading 
surface,  which  determines  the  loading  duration  by  the  first  wave  (precursor);  2)  duration  of 
the  loading  due  to  the  main  loading  wave;  and  3)  effect  of  the  precursor  rarefaction  from  the 
free  surface  of  the  sample. 

The  jump-like  approximation  is  used  for  argumentation  [47]  of  the  anomalous  behaviour  that 
includes: 

i)  assumption  of  a  much  larger  energy  increase  for  porous  material  (from  the  energy 
balance  equation  through  the  jump  in  system  (1))  due  to  a  larger  volume 
increment  when  compared  with  the  energy  increase  for  solid  material; 

ii)  use  of  Mie-Gruneisen  equation  of  state  for  both  a  solid  material  and  a  porous 
material  (mixture  of  gaseous  and  solid  phases),  using  a  single  Gruneisen 
coefficient;  and 

iii)  assumption  of  different  physical  mechanisms  of  the  material  response  for  porous 
and  solid  substances,  focusing  on  the  activation  of  a  high-temperature  mechanism 
in  a  porous  material  due  to  the  extreme  energy  increase;  the  only  obvious 
difference  between  these  substances  is  the  presence  of  the  gaseous  phase,  so  the 
high-temperature  mechanism  can  be  activated  only  by  compression  of  the  gaseous 
phase. 

When  the  jump-like  approximation  is  critically  reviewed  it  appears  that: 

i)  the  energy  increase  might  be  much  less  than  expected  due  to  several  shock 
transformations;  firstly,  pre-compacting  the  material  by  the  precursor  and, 
secondly,  transforming  it  into  the  final  pressure  state;  as  a  result,  the  energy 
calculated  from  system  (1)  is  incremented  in  several  steps  (in  a  more  isoentropic 
fashion); 

ii)  the  Mie-Gruneisen  approximation  can  be  valid  for  each  of  the  phases.  However, 
this  approximation  with  a  single  Gruneisen  coefficient  might  be  not  applicable  to 
a  mixture  of  the  two  phases,  where  the  gaseous  phase  can  change  this  coefficient 
significantly  during  the  shock  compression;  and 
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iii)  while  some  high-temperature  effects  are  obvious  for  porous  materials,  the  effect 
might  be  overestimated.  For  example,  data  [1,  5,  6]  claim  that  the  experimental 
data  for  the  porous  materials  in  the  atmospheric  air  environment  and  for  samples 
with  the  air  evacuated  from  the  pores  are  similar.  This  claim  does  actually  reduce 
the  high-temperature  contribution  from  compressed  air  in  the  pores,  at  least  into 
the  (D-D)-data. 

Commenting  on  the  first  item  above,  as  seen  from  data  [14-16],  it  is  hard  to  resolve  the 
compression  curve  structure  in  order  to  find  shock  wave  stability  conditions  similar  to  the 
phase  transition  or  elastic-plastic  transformations.  Nevertheless,  from  the  general  assessment 
of  the  precursor  velocity  it  is  seen  that  at  a  precursor  velocity  of  the  order  of  4  km/  s  for  the 
porous  copper  (that  is  not  in  contradiction  with  experiments  [1,6]),  wave  stability  could  be 
expected  at  the  behind-shock  pressures  of  more  than  26  GPa  and,  correspondingly,  at  the 
behind-shock  particle  velocities  of  more  than  3  km/  s.  This  is  in  agreement  with  the 
compression  behaviour  and  calculated  velocities.  It  is  remarkable,  that  the  abnormality  area  is 
within  the  possible  instability  area,  so  the  wave  splitting/  instability,  resulting  in  a  gradual 
energy  increments  that  are  lower  than  the  jump-like  energy  increase,  might  be  associated  with 
the  energy  overestimate. 


9.  Conclusions 


Direct  constitutive  modelling  has  been  used  to  simulate  Hugoniot  experiments  on  porous 
copper  with  a  porosity  m  =  4.  The  present  calculations  cannot  fully  reproduce  the  conditions 
of  the  experiments  because  the  published  experimental  set-ups  lack  the  necessary  details. 
Keeping  in  mind  significant  uncertainties  in  the  description  of  experimental  set-ups, 
parameters  of  the  set-ups,  and  the  procedures  of  data  IM-reduction  and  elimination  of  the 
shock  wave  attenuation  effects  applied  to  the  published  experimental  data,  the  calculations 
demonstrate  a  consistent  agreement  with  shock  velocity  data  obtained  in  the  experiments. 

The  present  work  shows  that  a  satisfactory  description  of  the  experiments  can  be  obtained 
without  assumptions  of  an  abnormal  behaviour  of  a  porous  material.  Resistance  of  material  to 
compression  due  to  possible  high-temperature  expansion  can  easy  be  overestimated  if  the 
constitutive  behaviour  of  porous  material  is  not  taken  into  account.  The  abnormality  may  also 
be  associated  with  the  wave  attenuation  effects  and  be  a  result  of  the  IM-reduction  procedure 
applied  to  the  data,  which  relates  the  Hugoniot  states  to  the  particle  velocity  of  material  at  the 
screen-sample  interface. 

Non-equilibrium  behind  the  shock  front  could  be  an  essential  factor  when  constructing  a 
Hugoniot  in  the  pressure-density  space.  The  Hugoniot  is  dependent  on  the  location  of  a  state 
on  the  shock  profile  where  the  state  is  taken  for  calculation  of  the  Hugoniot.  As  a  result,  for 
the  materials  exhibiting  distinct  constitutive  behaviour,  the  Hugoniot  may  depend  on  the 
sample  thickness  even  if  the  shock  velocity  does  not  change  essentially  with  the  thickness.  The 
present  work  has  attempted  to  evaluate  possible  contributions  of  the  constitutive  behaviour  to 
the  data  interpretation,  which  can  be  so  significant  as  to  appear  as  an  abnormality. 


30 


DSTO-TR-2026 


10.  References 


1.  Trunin  R.F.,  Shock  Compression  of  Condensed  Materials,  Cambridge  University  Press, 
Cambridge,  UK,  1998. 

2.  Walsh  J.M.,  Rice  M.H.,  McQueen  R.G.,  and  Yarger  F.L.,  Shock-Wave  Compressions  of 
Twenty-Seven  Metals.  Equations  of  State  of  Metals,  Phys.  Rev.,  108,  196-216  (1957), 
Issue  2-15  October  1957. 

3.  Al'tshuler  L.V.,  Krupnikov  K.K.,  and  Brazhnik  M.I.,  Dynamic  compressibility  of 
metals  under  pressures  from  400,000  to  4,000,000  atmospheres,  Soviet  Physics  JETP, 
1958,  v.  7,  n.  4,  pp.  614-619. 

4.  Al'tshuler  L.V.,  Development  of  Dynamical  Methods  of  Investigation  of  High 
Pressures  in  Russia,  In:  High-Pressure  Shock  Compression  of  Solids  VII,  edited  by  V  E 
Fortov,  L  V  Al'tshuler,  R  F  Trunin,  and  A I  Funtikov,  Springer-Verlag,  2004,  pp.  1-37. 

5.  Trunin  R.F.,  Krupnikov  K.K.,  Simakov  G.V.,  and  Funtikov  A.I.,  Shock-Wave 
Compression  of  Porous  Metals,  In:  High-Pressure  Shock  Compression  of  Solids  VII, 
edited  by  V  E  Fortov,  L  V  Al'tshuler,  R  F  Trunin,  and  A I  Funtikov,  Springer-Verlag, 
2004,  pp.  177-195. 

6.  Trunin  R.F.,  Simakov  G.V.,  Sutulov  Yu.N.,  Medvedev  A.B.,  Rogozkin  B.D.,  and 
Fedorov  Yu.E.,  Compressibility  of  porous  metals  in  shock  waves,  Soviet  Physics  JETP, 
1989,  v.  69,  n.  3,  pp.  580-588. 

7.  Trunin  R.F.,  Shock  compressibility  of  condensed  materials  in  strong  shock  waves 
generated  by  underground  nuclear  explosions,  Physics-Uspekhi,  1994,  v.  37,  n.  11,  pp. 
1123-1145. 

8.  Trunin  R.F.,  Medvedev  A.B.,  Funtikov  A.I.,  Podurets  M.A.,  Simakov  G.V.,  and 
Sevast'yanov  A.G.,  Shock  compression  of  porous  iron,  copper,  and  tungsten,  and  their 
equation  of  state  in  the  terapascal  pressure  range,  Soviet  Physics  JETP,  1989,  v.  68,  n.  2, 
pp.  356-361. 

9.  Kormer  S.B.,  Funtikov  A.I.,  Urlin  V.D.,  and  Kolesnikova  A.N.,  Dynamic  compression 
of  porous  metals  and  the  equation  of  state  with  variable  specific  heat  at  high 
temperatures,  Soviet  Physics  JETP,  1962,  v.  15,  n.  3,  pp.  477-488. 

10.  Zubarev  V.N.,  Podurets  M.A.,  Popov  L.V.,  Simakov  G.V.,  and  Trunin  R.F.,  Shock 
compressibility  and  equation  of  state  for  copper  in  a  high  pressure  region  [in  Russian], 
In:  Detonation:  Critical  phenomena.  Physical-Chemical  transformations  in  shock 
waves  (F.I.  Dubovitsky,  ed.).  Dept,  of  Institute  of  Chemical  Physics,  USSR  Academy  of 
Sciences,  Chernogolovka,  1978,  pp.  61-65. 

11.  Marsh  S.P.  (Ed.),  LASL  Shock  Hugoniot  Data,  University  of  California  Press,  Berkley 
and  Los  Angeles,  California,  1980. 

12.  van  Thiel  M.,  Shaner  J.,  and  Salinas  E.  (Eds.),  Compendium  of  Shock  Wave  Data, 
Lawrence  Livermore  Laboratory  Report  UCRL-50108,  University  of  California, 
Livermore,  Vol.  1, 1977. 

13.  Bolkhovitinov  L.C.  and  Khvostov  Yu.B.,  Polystyrene  and  copper  isentropes  obtained 
on  the  data  of  shock  compressed  materials  of  high  porosity,  Lett.  Appl.  Engng.  Sci., 
1984,  v.  22,  n.  4,  pp.  491-504. 

14.  Boade  R.R.,  Compression  of  Porous  Copper  by  Shock  Waves,  J.  Appl.  Physics,  1968,  v. 
39,  n.  12,  pp.  5693-5702. 


31 


DSTO-TR-2026 


15.  Boade  R.R.,  Principle  Hugoniot,  Second-Shock  Hugoniot,  and  Release  Behavior,  of 
Pressed  Copper  Powder,  J.  Appl.  Physics,  1970,  v.  41,  n.  11,  pp.  4542-4551. 

16.  Linde  R.K.,  Seaman  L.,  and  Schmidt  D.N.,  Shock  response  of  porous  copper,  iron, 
tungsten,  and  polyurethane,  J.  Appl.  Physics,  1972,  v.  43,  n.  8,  pp.  3367-3375. 

17.  Trunin  R.F.,  Podurets  M.A.,  Moiseev  B.N.,  Simakov  G.V.,  and  Popov  L.V.,  Relative 
compressibility  of  copper,  cadmium  and  lead  at  high  pressures,  Soviet  Physics  JETP, 
1969,  v.  29,  n.  4,  pp.  630-631. 

18.  Alekseev  Yu.L.,  Ratnikov  V.P.,  and  Rybakov  A.P.,  Shock  adiabats  of  porous  metals,  J. 
Appl.  Mech.  Techn.  Physics,  1971,  v.  12,  n.  2,  pp.  257-262. 

19.  Bakanova  A.  A.,  Dudoladov  I.P.,  and  Sutulov  Yu.N.,  Shock  compressibility  of  porous 
tungsten,  molybdenum,  copper,  and  aluminium  in  the  low  pressure  domain,  J.  Appl. 
Mech.  Techn.  Physics,  1974,  v.  15,  n.  2,  pp.  241-245. 

20.  Al'tshuler  L.V.,  Bushman  A.V.,  Zhernokletov  M.V.,  Zubarev  V.N.,  Leont'ev  A.  A.,  and 
Fortov  V.E.,  Unloading  isentropes  and  the  equation  of  state  of  metals  at  high  energy 
densities,  Soviet  Physics  JETP,  1980,  v.  51,  n.  2,  pp.  373-383. 

21.  Gryaznov  V.K.,  Fortov  V.E.,  Zhernokletov  M.V.,  Simakov  G.V.,  Trunin  R.F.,  Trusov 
L.I.,  and  Iosilevskii  I.L.,  Shock  compression  and  thermodynamics  of  highly  nonideal 
metallic  plasma,  Soviet  Physics  JETP,  1998,  v.  87,  n.  4,  pp.  678-690. 

22.  Nellis  W.J.,  Moriarty  J.A.,  Mitchell  A.C.,  Ross  M.,  Dandrea  R.G.,  Ashcroft  N.W., 
Holmes  N.C.,  and  Gathers  G.R.,  Metals  Physics  at  Ultrahigh  Pressure:  Aluminum, 
Copper,  and  Lead  as  Prototypes,  Phys.  Rev.  Lett.,  1988,  v.  60,  n.  14,  pp.  1414-1417. 

23.  Butcher  B.  M.,  Carroll  M.  M.,  and  Holt  A.C.,  Shock-wave  compaction  of  porous 
aluminum,  J.  Appl.  Phys.,  v.  45,  n.  9, 1974,  pp.  3864-3875. 

24.  Resnyansky  A.D.  and  Bourne  N.K.,  Shock-wave  compression  of  a  porous  material,  J 
Appl.  Physics,  2004,  v.  95,  n.  4,  pp.  1760-1769. 

25.  Trunin  R.F.,  Shock  compression  of  condensed  materials  (laboratory  studies),  Physics- 
Uspekhi,  2001,  v.  44,  n.  4,  pp.  371-396. 

26.  Al'tshuler  L.V.,  Trunin  R.F.,  Krupnikov  K.K.,  and  Panov  N.V.,  Explosive  laboratory 
devices  for  shock  wave  compression  studies,  Physics-Uspekhi,  1996,  v.  39,  n.  5,  pp. 
539-544. 

27.  Al'tshuler  L.V.,  Trunin  R.F.,  Urlin  V.D. ,  Fortov  V.E.,  and  Funtikov  A.I.,  Development 
of  dynamic  high-pressure  techniques  in  Russia,  Physics-Uspekhi,  1999,  v.  42,  n.  3,  pp. 
261-280. 

28.  Asay  J.R.,  Shock  wave  paradigms  and  new  challenges.  In:  Shock  Compression  of 
Condensed  Matter  -  2001,  edited  by  M  D  Furnish,  N  N  Thadhani,  and  Y  Horie,  (AIP, 
New  York,  2002),  CP620,  2002,  pp.  26-35. 

29.  Mallory  H.D.,  Propagation  of  Shock  Waves  in  Aluminum,  J  Appl.  Physics,  1955,  v.  26, 
n.  5,  pp.  555-559. 

30.  Goranson  R.W.,  Bancroft  D.,  Burton  B.L.,  Blechar  T.,  Houston  E.E.,  Gittings  E.F.,  and 
Landeen  S.A.,  Dynamic  Determination  of  the  Compressibility  of  Metals,  J  Appl. 
Physics,  1955,  v.  26,  n.  12,  pp.  1472-1479. 

31.  Bugaeva  V.A.,  Evstigneev  A.A.,  and  Trunin  R.F.,  Analysis  of  Calculation  Data  on  the 
Adiabats  of  Expansion  for  Copper,  Iron,  and  Aluminium,  High  Temperature,  1996,  v. 
34,  n.  5,  pp.  674-680. 

32.  Avrorin  E.N.,  Vodolaga  B.K.,  Simonenko  V.A.,  and  Fortov  V.E.,  Intense  shock  waves 
and  extreme  states  of  matter,  Physics-Uspekhi,  1993,  v.  36,  n.  5,  pp.  337-364. 


32 


DSTO-TR-2026 


33.  Grady  D.E.,  Hydrodynamic  compressibility  of  silicon  carbide  through  shock 
compression  of  metal-ceramic  mixtures,  J  Appl.  Physics,  1994,  v.  75,  n.  1,  pp.  197-202. 

34.  Grady  D.,  Analytical  solutions  and  constitutive  relations  for  shock  propagation  in 
porous  media.  In:  Shock  Compression  of  Condensed  Matter  -  2003,  edited  by  M  D 
Furnish,  Y  M  Gupta,  and  J  W  Forbes,  (AIP,  New  York,  2004),  CP706, 2004,  pp.  205-208. 

35.  Resnyansky  A.D.  and  Bourne  N.K.,  Shock  compression  of  dry  and  hydrated  sand.  In: 
Shock  Compression  of  Condensed  Matter  -  2003,  edited  by  M  D  Furnish,  Y  M  Gupta, 
and  J  W  Forbes,  (AIP,  New  York,  2004),  CP706,  2004,  pp.  1474-1477. 

36.  Carroll  M.M.  and  Holt  A.C.,  Static  and  Dynamic  Pore-Collapse  Relations  for  Ductile 
Porous  Materials,  J  Appl.  Physics,  1972,  v.  43,  n.  4,  pp.  1626-1636. 

37.  Kapila  A.K.,  Menikoff  R.,  Bdzil  J.B.,  Son  S.F.,  and  Stewart  D.S.,  Two-phase  modeling  of 
deflagration-to-detonation  transition  in  granular  materials:  Reduced  equations. 
Physics  of  Fluids,  2001,  v.  13,  n.  10,  pp.  3002-3024. 

38.  Godunov  S.K.,  Kozin  N.S.,  and  Romenskii  E.I.,  Equations  of  state  of  the  elastic  energy 
of  metals  in  the  case  of  a  nonspherical  strain  tensor,  J.  Appl.  Mech.  and  Techn.  Phys., 
1974,  v.  15,  n.  2,  pp.  246-250. 

39.  Resnyansky  A.D.  and  Romensky  E.I.,  Numerical  Modelling  the  Non-Equilibrium 
Phase  Transitions  in  the  Continuum  Media  By  a  Nonlinear  Kinetics,  DYMAT  Journal, 
1995,  v.  2,  n.  3/4,  pp.  223-228. 

40.  Godunov  S.K.,  Zabrodin  A.V.,  Ivanov  M.Ya.,  Kraiko  A.N.,  and  Prokopov  G.P., 
Numerical  Solution  of  Multi-Dimensional  Problems  of  Gas  Dynamics,  Nauka  Press, 
Moscow,  1976  (French  transl:  Resolution  Numerique  des  Problemes 
Multidimensionnels  de  la  Dynamique  des  Gaz,  Mir,  Moscow,  1979. 

41.  Resnyansky  A.D.,  A  technique  of  use  of  arbitrary  equations  of  state  in  the  Riemann 
problem  for  the  Godunov  scheme  [in  Russian],  Dynamics  of  Continua  (Dinamika 
sploshnoi  sredy),  Favrentyev  Institute  of  Hydrodynamics,  Novosibirsk,  USSR,  1988,  n. 
86,  pp.  88-95. 

42.  Godunov  S.K.  and  Romenskii  E.I.,  "Elements  of  Continuum  Mechanics  and 
Conservation  Faws",  Kluwer  Academic  Publ.,  N.Y.,  2003. 

43.  Resnyansky  A.D.,  DYNA-modelling  of  the  high-velocity  impact  problems  with  a  split- 
element  algorithm,  Int.  J.  Impact  Eng.,  2002,  v.  27,  n.  7,  pp.  709-727. 

44.  Trunin  R.F.  and  Panov  N.V.,  Shock  Compression  of  Porous  Copper  at  Megabar 
Pressures,  High  Temperature,  2000,  v.  38,  n.  5,  pp.  754-758. 

45.  Findholm  U.S.  and  Yeakley  F.M.,  High  Strain-rate  Testing:  Tension  and  Compression, 
1968,  Exp.  Mech.,  v.  8,  pp.  1-9. 

46.  Heyda  J.F.,  Plate-gap  model  of  a  porous  solid  and  its  application  to  impact  by  reduced 
density  projectiles,  NASA  report  NASA  CR-1140,  NASA-Fangley,  1968. 

47.  Zel'dovich  Ya.  B.  and  Raizer  Yu.  P.,  "Physics  of  shock  waves  and  high-temperature 
hydrodynamic  phenomena",  vol.  2,  Academic  Press,  N.Y.,  1967. 

48.  Godunov  S.K.,  Zabrodin  A. V.,  Ivanov  M.Y.,  Kraiko  A.N.,  Prokopov  G.P.,  "Resolution 
numerique  des  problemes  multidimensionnels  de  la  dynamique  des  gaz",  Mir, 
Moscow,  1979  [in  French,  transl.  from  "Numerical  solution  of  multi-dimensional 
problems  in  gas  dynamics",  Nauka  Press,  Moscow,  1976,  in  Russian]. 

49.  Kolgan  V.P.,  Application  of  the  minimum-derivative  principle  in  the  construction  of 
finite-difference  schemes  for  numerical  analysis  of  discontinuous  solutions  in  gas 
dynamics  [in  Russian],  Uch.  Zap.  TsaGI,  1972,  v.  3,  n.  6,  pp.  68-77. 


33 


DSTO-TR-2026 


34 


DSTO-TR-2026 


Appendix  A:  Nonlinear  Riemann  problem  for  arbitrary 

equation  of  state 


This  Appendix  describes  briefly  the  method  of  calculation  of  the  Riemann  problem  for  an 
isotropic  material  with  an  arbitrary  two-parameter  Equation  Of  State  (EOS);  the  method  is 
applicable  to  the  Godunov  scheme  with  the  second  order  spatial  accuracy.  This  method  has 
been  used  for  a  variety  of  two-phase  materials  and  the  mathematical  rationale  of  the  method 
has  been  published  earlier  in  [41] .  The  method  allows  one  to  approximate  an  arbitrary  EOS  as 
an  extension  of  the  ideal  gas  EOS  in  the  following  form: 


p  =  A(S 


V 

P_ 

Pi)  j 


(19) 


we  require  here  fulfilment  of  the  condition  A'  >  0.  Suppose  we  have  an  arbitrary  EOS  in  the 
following  form 

P  =  f(p,S)  ,  (20) 


where  the  function/obeys  the  following  conditions  of  normal  material:/p>0  and/pp>0.  For  the 
Riemann  problem  with  initial  parameters  pi,  pi,  Ui  and  pi,  p2,  Ui  on  the  left  and  right  from  a 
contact  surface  of  a  'shock  tube'  (we  denote  the  materials  with  such  the  initial  states  as 
materials  '1'  and  '2')  we  need  to  find  a  solution  at  the  contact  surface  at  t  >  0.  It  should  be 
noted  that  from  the  EOS  (20)  the  parameters  Si  and  S2  at  t  =  0  can  be  found  as  well.  Exact 
solutions  of  the  Riemann  problem  for  the  Godunov  scheme  are  documented  well  [48]  for  the 
EOS  in  the  form  (19)  even  if  the  constants  po  and  y  are  different  in  the  materials  '1'  and  '2'.  In 
the  present  work,  the  Godunov  scheme  [48]  is  used  in  the  second  order  accuracy  modification 
published  in  [49]. 


Now  for  the  materials  '1'  and  '2'  we  can  introduce  the  following  constants  associated  with 
the  current  states  (pi,  pi.  Si)  and  (p2,  pi,  Si): 


Yi  =1  + 

Yi  =1  + 


PiPppiPiA  1) 
pp{pi>s  1) 

PiPppiPi’Si) 

Pp{p2’S2) 


0  PlPp{Pl’S  l)  ~Y 

Px  = - -P\Pi,Sx), 

n 

0  P2P piPi’Si)  — /  c  \ 
Pi  = - P\Pi,S2). 

Yi 


Then,  we  can  use  the  EOS  (19)  with  the  constants  yi  and  pp  instead  of  y  and  p°  for  the  material 
'1'  and  the  same  EOS  with  the  constants  \i  and  pi0  for  the  material  '2'.  It  has  been  proved  in 
the  publication  [41]  that  solution  of  the  Riemann  problem  using  these  EOSs  approximates 
with  second  order  accuracy  the  solution  of  the  Riemann  problem  using  the  EOS  (20). 
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Appendix  B:  Impedance  mismatch  procedure 

This  Appendix  describes  the  procedure  of  calculation  of  the  shock  Hugoniot  parameters  using 
the  Impedance  Mismatch  (IM)  method  and  evaluates  possible  errors  for  the  parameters  due  to 
errors  in  the  shock  velocity  determination.  Assuming  the  jump-like  approximation  for  the 
shock  wave  and  using  the  standard  assumptions  for  a  Hugoniot  [3],  we  relate  the  shock 
velocity  to  the  particle  velocities  in  the  standard  material  via  a  linear  dependence: 


Ds=  a  +  b  U  .  (21) 

Referring  to  the  schematic  in  Fig.  2(b),  the  Hugoniot  Hs  is  approximated  as  the  following 
quadratic  dependence 

P  =  p0U  \a  +  b  U)  ,  (22) 

when  neglecting  the  initial  pressure  and  particle  velocity  in  front  of  the  shock.  Then,  in 
agreement  with  [3],  we  approximate  the  rarefaction  curve  Rs  by  the  mirrored  Hugoniot  Hs 
with  respect  to  the  impact  velocity  2Uo  that  is  expressed  by  the  following  equation 

P  =  -  pf  (U  -  2U0)  ip  -b\ U  -  2U0))  .  (23) 

It  should  be  noted  that  the  particle  velocity  in  the  screen  Uo  can  be  found  from  the  Hugoniot 
relation  (21)  if  the  shock  velocity  Ds  in  the  standard  has  been  recorded  in  the  experiment. 
When  the  shock  velocity  D  (Dp  in  Fig.  2)  in  the  investigated  material  has  also  been  measured, 
the  particle  velocity  Uh  at  point  2  in  Fig.  2(b)  can  be  calculated  from  the  intersection  of  the 
Hugoniot  Hs  defined  by  equation  (23)  with  the  Rayleigh  line  P  =  poo  'U  D.  This  gives 


Uh  =  2  U0 


+ 


Psa  +  PoaD  -  ^{psa  +  P00DY  +Sp.M 


00DU0 


l(2psb). 


(24) 


Suppose  we  have  tabulated  the  (D-U)-Hugoniot  (21)  of  the  standard  with  very  good  accuracy 
so  that  we  can  neglect  inaccuracies  in  determination  of  the  curve.  Thus,  we  can  also  neglect 
inaccuracies  in  Uothat  is  calculated  from  (21)  and  appears  in  (24). 

Let  us  evaluate  the  possible  error  for  the  velocity  Uh  if  an  100x5%  error  has  been  made  in  the 
determination  of  the  shock  velocity  D.  Denote  the  erroneous  shock  velocity  by  D'  so  that 

D'-D  =  5-D  (25) 

Then,  the  corresponding  particle  velocity  Uh'  is  obtained  from  (24)  with  D  replaced  by  D'. 
Subtracting  Uh'  from  Uh  (obtained  from  (24)),  the  error  estimate  for  Uh  can  be  found  as  follows 

U’h-Uh  D'-D  ^  ^(Pfa  +  PMD)2  +8pfbp00DU0 

- = - ,  A  = - .  (26) 

Uh  KD  PoqD 
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Assuming  that  the  Hugoniot  Hs  is  built  up  without  errors  and  the  locus  of  the  rarefaction 
states  in  Fig.  2  follows  the  curve  Hs  ,  it  is  seen  that  K  >  1  and  the  error  in  Uh  is  of  the  same 
order  as  the  error  in  D.  Thus,  from  (25)  and  (26)  an  error  for  the  disturbed  particle  velocity  is 
Uh'  -  Uh  =  5  Uh/ K.  Thus,  the  error  in  the  particle  velocity  is  K  times  less  than  the  error  in  the 
shock  velocity,  and  the  K-coefficient  depends  mainly  on  the  initial  parameters,  the  Hugoniot 
of  the  standard  material  and  the  magnitude  of  the  wave  propagating  in  the  standard  and  test 
materials  (Uo  and  D).  For  the  aluminium  standard  and  the  porous  material  analysed  this 
coefficient  varies  from  10  at  a  low  velocity  of  impact  down  to  3  at  the  highest  impact  velocity 
from  the  load  range  analysed.  It  should  be  noted  that  this  is  an  upper  estimate  of  the  error 
because  at  the  high  velocity  of  impact  the  mirror  approximation  of  Rs  by  Hs  deteriorates  and 
at  the  low  velocity  of  impact  the  precursor  multiple  reflection  and  circulation  effects  interfere 
with  the  approximation.  However,  the  Uh-error  is  at  least  not  larger  than  the  D-error.  The 
pressure  error,  in  the  first  order,  is 

P'  -  P  =  poo(H'Uh'  -  DUh)  “  8  -P  '(1+1/ K)  .  (27) 


Thus,  the  P-error  can  be  slightly  greater  than  the  D-error  although  of  similar  magnitude.  The 
last  parameter  we  estimate  is  error  in  determination  of  density.  Evaluating  density  from  poo/  p 
=  (D  -  Uh)/D,  we  have 


P oo  _  P  oo 
P’  P 


u,  u, 

5 

D  D’ 


K. 
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1  + 


n 

K , 


(28) 


For  the  porous  copper  D  ~  1.5Uh  and,  thus,  the  coefficient  Kr  is  of  the  order  of  1  (more 
precisely,  between  0.75  and  1,  neglecting  the  non-equilibrium  and  the  mirror  Hugoniot 
issues).  When  calculating  density,  we  have 


P~P’  P 

P  Poo 


■  K,.S . 


(29) 


This  means  that  at  compression  of  a  porous  material  when  the  density  approaches  the  density 
of  the  solid  phase  po,  the  p-error  may  be  as  high  as  mKr5  (m  =  po/ poo).  For  the  porous  copper 
considered  above  m  =  4  and  the  error  varies  between  35  and  45,  thus,  exceeding  the  D-error 
by  a  factor  of  3-4. 
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